!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                                                             
!   glide_types.F90 - part of the Community Ice Sheet Model (CISM)  
!                                                              
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!   Copyright (C) 2005-2018
!   CISM contributors - see AUTHORS file for list of contributors
!
!   This file is part of CISM.
!
!   CISM is free software: you can redistribute it and/or modify it
!   under the terms of the Lesser GNU General Public License as published
!   by the Free Software Foundation, either version 3 of the License, or
!   (at your option) any later version.
!
!   CISM is distributed in the hope that it will be useful,
!   but WITHOUT ANY WARRANTY; without even the implied warranty of
!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!   Lesser GNU General Public License for more details.
!
!   You should have received a copy of the Lesser GNU General Public License
!   along with CISM. If not, see <http://www.gnu.org/licenses/>.
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif

module glide_types

  !> Holds type definitions for the derived types used by each 
  !> instance of the ice model. Originally, each of these types
  !> was a module containing variables, which were used as containers
  !> for global variables. However, the need to allow for multiple
  !> ice model instances meant that the nested derived types were instituted
  !> instead. However, there is probably one too many levels in this scheme. 
  !> It would be better if the different types here were contained in the 
  !> higher-level instance type (\texttt{glint\_instance}), rather than 
  !> the intermediate model type (\texttt{glide\_global\_type}). 
  !> 
  !> Note that this \emph{is} now where the defaults are defined for these
  !> variables.

!TODO - Clean up the glide_global type so it holds fewer subtypes?
!       For example, we could replace some work types (tempwk, velowk) with local arrays and parameters.

  use glimmer_sparse_type
  use glimmer_global, only: sp, dp
  use glimmer_ncdf
  use profile
  use glimmer_coordinates, only: coordsystem_type
  use glimmer_map_types
  use glimmer_physcon
  use glimmer_paramets, only: unphys_val

  implicit none

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! Constants that describe the options available
  ! We use these integer parameters elsewhere in the code to avoid
  !  hardwiring of option numbers

  ! basic Glimmer/Glide options

  integer, parameter :: GLOBAL_BC_PERIODIC = 0  ! doubly periodic
  integer, parameter :: GLOBAL_BC_OUTFLOW = 1   ! free outflow; scalars in global halo set to zero
  integer, parameter :: GLOBAL_BC_NO_PENETRATION = 2   ! no penetration; outflow set to zero at global boundaries
                                                       ! NOTE: The no-penetration option is currently supported for Glissade only

  integer, parameter :: DYCORE_GLIDE = 0     ! old shallow-ice dycore from Glimmer
!WHL - The Glam dycore is no longer supported; choosing whichdycore = 1 will give a fatal error.
!      Might repurpose the '1' slot for a future version of Glissade.
  integer, parameter :: DYCORE_GLAM = 1      ! Payne-Price finite-difference solver
  integer, parameter :: DYCORE_GLISSADE = 2  ! prototype finite-element solver
  integer, parameter :: DYCORE_ALBANYFELIX = 3  ! External Albany-Felix finite-element solver
  integer, parameter :: DYCORE_BISICLES = 4     ! BISICLES-Chombo external FVM solver

  integer, parameter :: EVOL_PSEUDO_DIFF = 0    ! glide only
  integer, parameter :: EVOL_ADI = 1            ! glide only
  integer, parameter :: EVOL_DIFFUSION = 2      ! glide only
  integer, parameter :: EVOL_INC_REMAP = 3      ! glissade only
  integer, parameter :: EVOL_UPWIND = 4         ! glissade only
  integer, parameter :: EVOL_NO_THICKNESS = 5   ! glissade only

  !NOTE: Use option 1 for prognostic temperature with any dycore
  !      Option 3 is under construction

  integer, parameter :: TEMP_SURFACE_AIR_TEMP = 0
  integer, parameter :: TEMP_PROGNOSTIC = 1
  integer, parameter :: TEMP_STEADY = 2
  integer, parameter :: TEMP_ENTHALPY = 3

  integer, parameter :: TEMP_INIT_ZERO = 0
  integer, parameter :: TEMP_INIT_ARTM = 1
  integer, parameter :: TEMP_INIT_LINEAR = 2
  integer, parameter :: TEMP_INIT_ADVECTIVE_DIFFUSIVE = 3
  integer, parameter :: TEMP_INIT_EXTERNAL = 4

  integer, parameter :: FLWA_CONST_FLWA = 0
  integer, parameter :: FLWA_PATERSON_BUDD_CONST_TEMP = 1
  integer, parameter :: FLWA_PATERSON_BUDD = 2
  integer, parameter :: FLWA_INPUT = 3

  integer, parameter :: BTRC_ZERO = 0
  integer, parameter :: BTRC_CONSTANT = 1
  integer, parameter :: BTRC_CONSTANT_BWAT = 2
  integer, parameter :: BTRC_CONSTANT_BPMP = 3
  integer, parameter :: BTRC_LINEAR_BMLT = 4
  integer, parameter :: BTRC_TANH_BWAT = 5

  integer, parameter :: BWATER_NONE  = 0
  integer, parameter :: BWATER_LOCAL = 1
  integer, parameter :: BWATER_FLUX  = 2
  integer, parameter :: BWATER_CONST = 3
  ! option 4 is deprecated; if selected, the code will throw a fatal error
  integer, parameter :: BWATER_OCEAN_PENETRATION = 4

  integer, parameter :: BMLT_FLOAT_NONE = 0
  integer, parameter :: BMLT_FLOAT_MISMIP = 1
  integer, parameter :: BMLT_FLOAT_CONSTANT = 2
  integer, parameter :: BMLT_FLOAT_DEPTH = 3
  integer, parameter :: BMLT_FLOAT_EXTERNAL = 4
  integer, parameter :: BMLT_FLOAT_MISOMIP = 5

  integer, parameter :: BASAL_MBAL_NO_CONTINUITY = 0
  integer, parameter :: BASAL_MBAL_CONTINUITY = 1

  integer, parameter :: SMB_INPUT_MYR_ICE = 0     ! use 'acab' for input
  integer, parameter :: SMB_INPUT_MMYR_WE = 1     ! use 'smb' for input

  integer, parameter :: OVERWRITE_ACAB_NONE = 0
  integer, parameter :: OVERWRITE_ACAB_ZERO_ACAB = 1
  integer, parameter :: OVERWRITE_ACAB_THCKMIN = 2

  integer, parameter :: GTHF_UNIFORM = 0
  integer, parameter :: GTHF_PRESCRIBED_2D = 1
  integer, parameter :: GTHF_COMPUTE = 2

  integer, parameter :: RELAXED_TOPO_DEFAULT = 0  ! topo and relx are separate input fields
  integer, parameter :: RELAXED_TOPO_INPUT = 1    ! set relx to input topg
  integer, parameter :: RELAXED_TOPO_COMPUTE = 2  ! Input topo in isostatic equilibrium
                                                  ! compute relaxed topo

  integer, parameter :: ISOSTASY_NONE = 0
  integer, parameter :: ISOSTASY_COMPUTE = 1

  integer, parameter :: LITHOSPHERE_LOCAL = 0
  integer, parameter :: LITHOSPHERE_ELASTIC = 1

  integer, parameter :: ASTHENOSPHERE_FLUID = 0
  integer, parameter :: ASTHENOSPHERE_RELAXING = 1

  integer, parameter :: CALVING_NONE = 0
  integer, parameter :: CALVING_FLOAT_ZERO = 1
  integer, parameter :: CALVING_FLOAT_FRACTION = 2
  integer, parameter :: CALVING_RELX_THRESHOLD = 3
  integer, parameter :: CALVING_TOPG_THRESHOLD = 4
  integer, parameter :: CALVING_GRID_MASK = 5
  integer, parameter :: CALVING_THCK_THRESHOLD = 6
  integer, parameter :: EIGENCALVING = 7
  integer, parameter :: CALVING_DAMAGE = 8
  integer, parameter :: CALVING_HUYBRECHTS = 9

  integer, parameter :: CALVING_INIT_OFF = 0
  integer, parameter :: CALVING_INIT_ON = 1

  integer, parameter :: CALVING_DOMAIN_OCEAN_EDGE = 0
  integer, parameter :: CALVING_DOMAIN_EVERYWHERE = 1

  integer, parameter :: VERTINT_STANDARD = 0
  integer, parameter :: VERTINT_KINEMATIC_BC = 1

  integer, parameter :: DM_DT_DIAG_KG_S = 0
  integer, parameter :: DM_DT_DIAG_GT_Y = 1

  integer, parameter :: DIAG_MINTHCK_ZERO = 0
  integer, parameter :: DIAG_MINTHCK_THKLIM = 1

  integer, parameter :: SIGMA_COMPUTE_GLIDE = 0
  integer, parameter :: SIGMA_EXTERNAL = 1
  integer, parameter :: SIGMA_CONFIG = 2
  integer, parameter :: SIGMA_COMPUTE_EVEN = 3
  integer, parameter :: SIGMA_COMPUTE_PATTYN = 4

  integer, parameter :: RESTART_FALSE = 0
  integer, parameter :: RESTART_TRUE = 1

  integer, parameter :: RESTART_EXTEND_VELO_FALSE = 0
  integer, parameter :: RESTART_EXTEND_VELO_TRUE = 1
  
  !basal proc option disabled for now
  integer, parameter :: BAS_PROC_DISABLED = 0
!!  integer, parameter :: BAS_PROC_FULLCALC = 1
!!  integer, parameter :: BAS_PROC_FASTCALC = 2

  ! higher-order options

  integer, parameter :: HO_EFVS_CONSTANT = 0
  integer, parameter :: HO_EFVS_FLOWFACT = 1
  integer, parameter :: HO_EFVS_NONLINEAR = 2

  integer, parameter :: HO_DISP_NONE = -1
  integer, parameter :: HO_DISP_SIA = 0
  integer, parameter :: HO_DISP_FIRSTORDER = 1

  integer, parameter :: HO_THERMAL_BEFORE_TRANSPORT = 0
  integer, parameter :: HO_THERMAL_AFTER_TRANSPORT = 1
  integer, parameter :: HO_THERMAL_SPLIT_TIMESTEP = 2
  
  !TODO - Deprecate the last two options? Rarely if ever used.
  integer, parameter :: HO_BABC_BETA_CONSTANT = 0
  integer, parameter :: HO_BABC_BETA_BPMP = 1
  integer, parameter :: HO_BABC_YIELD_PICARD = 2
  integer, parameter :: HO_BABC_PSEUDO_PLASTIC = 3
  integer, parameter :: HO_BABC_BETA_LARGE = 4
  integer, parameter :: HO_BABC_BETA_EXTERNAL = 5
  integer, parameter :: HO_BABC_NO_SLIP = 6
  integer, parameter :: HO_BABC_YIELD_NEWTON = 7
  integer, parameter :: HO_BABC_ISHOMC = 8
  integer, parameter :: HO_BABC_POWERLAW = 9
  integer, parameter :: HO_BABC_COULOMB_FRICTION = 10
  integer, parameter :: HO_BABC_COULOMB_POWERLAW_SCHOOF = 11
  integer, parameter :: HO_BABC_COULOMB_POWERLAW_TSAI = 12
  integer, parameter :: HO_BABC_POWERLAW_EFFECPRESS = 13
  integer, parameter :: HO_BABC_SIMPLE = 14

  integer, parameter :: HO_INVERSION_NONE = 0
  integer, parameter :: HO_INVERSION_COMPUTE = 1
  integer, parameter :: HO_INVERSION_PRESCRIBE = 2

  integer, parameter :: HO_BWAT_NONE = 0
  integer, parameter :: HO_BWAT_CONSTANT = 1
  integer, parameter :: HO_BWAT_LOCAL_TILL = 2

  !TODO - Remove option 2? Rarely used
  integer, parameter :: HO_EFFECPRESS_OVERBURDEN = 0
  integer, parameter :: HO_EFFECPRESS_BPMP = 1
  integer, parameter :: HO_EFFECPRESS_BMLT = 2
  integer, parameter :: HO_EFFECPRESS_OCEAN_PENETRATION = 3
  integer, parameter :: HO_EFFECPRESS_BWAT = 4

  integer, parameter :: HO_NONLIN_PICARD = 0
  integer, parameter :: HO_NONLIN_JFNK = 1

  integer, parameter :: HO_RESID_MAXU = 0
  integer, parameter :: HO_RESID_MAXU_NO_UBAS = 1
  integer, parameter :: HO_RESID_MEANU = 2
  integer, parameter :: HO_RESID_L2NORM = 3
  integer, parameter :: HO_RESID_L2NORM_RELATIVE = 4

  integer, parameter :: HO_SPARSE_PCG_INCH = -1
  integer, parameter :: HO_SPARSE_BICG = 0
  integer, parameter :: HO_SPARSE_GMRES = 1
  integer, parameter :: HO_SPARSE_PCG_STANDARD = 2
  integer, parameter :: HO_SPARSE_PCG_CHRONGEAR = 3
  integer, parameter :: HO_SPARSE_TRILINOS = 4

  integer, parameter :: HO_APPROX_LOCAL_SIA = -1
  integer, parameter :: HO_APPROX_SIA = 0
  integer, parameter :: HO_APPROX_SSA = 1
  integer, parameter :: HO_APPROX_BP = 2
  integer, parameter :: HO_APPROX_L1L2 = 3
  integer, parameter :: HO_APPROX_DIVA = 4

  integer, parameter :: HO_PRECOND_NONE = 0
  integer, parameter :: HO_PRECOND_DIAG = 1
  integer, parameter :: HO_PRECOND_SIA  = 2

  integer, parameter :: HO_GRADIENT_CENTERED = 0
  integer, parameter :: HO_GRADIENT_UPSTREAM1 = 1
  integer, parameter :: HO_GRADIENT_UPSTREAM2 = 2

  integer, parameter :: HO_GRADIENT_MARGIN_LAND = 0
  integer, parameter :: HO_GRADIENT_MARGIN_HYBRID = 1
  integer, parameter :: HO_GRADIENT_MARGIN_MARINE = 2

  integer, parameter :: HO_VERTICAL_REMAP_FIRST_ORDER = 0
  integer, parameter :: HO_VERTICAL_REMAP_SECOND_ORDER = 1

  integer, parameter :: HO_ASSEMBLE_BETA_STANDARD = 0
  integer, parameter :: HO_ASSEMBLE_BETA_LOCAL = 1

  integer, parameter :: HO_ASSEMBLE_TAUD_STANDARD = 0
  integer, parameter :: HO_ASSEMBLE_TAUD_LOCAL = 1

  integer, parameter :: HO_ASSEMBLE_BFRIC_STANDARD = 0
  integer, parameter :: HO_ASSEMBLE_BFRIC_LOCAL = 1

  integer, parameter :: HO_CALVING_FRONT_NO_SUBGRID = 0
  integer, parameter :: HO_CALVING_FRONT_SUBGRID = 1

  integer, parameter :: HO_GROUND_NO_GLP = 0
  integer, parameter :: HO_GROUND_GLP = 1
  integer, parameter :: HO_GROUND_GLP_QUADRANTS = 2

  integer, parameter :: HO_GROUND_BMLT_NO_GLP = 0
  integer, parameter :: HO_GROUND_BMLT_FLOATING_FRAC = 1
  integer, parameter :: HO_GROUND_BMLT_ZERO_GROUNDED = 2

  integer, parameter :: HO_FLOTATION_FUNCTION_PATTYN = 0
  integer, parameter :: HO_FLOTATION_FUNCTION_INVERSE_PATTYN = 1
  integer, parameter :: HO_FLOTATION_FUNCTION_LINEAR = 2

  integer, parameter :: HO_ICE_AGE_NONE = 0 
  integer, parameter :: HO_ICE_AGE_COMPUTE = 1 

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_general

    !> Holds fundamental parameters of the ice model geometry.

    integer :: ewn = 0  !> The number of grid-points in the E-W direction.
    integer :: nsn = 0  !> The number of grid-points in the N-S direction.
    integer :: upn = 3  !> The number of vertical levels in the model.

    type(coordsystem_type) :: ice_grid  !> coordinate system of the ice grid
    type(coordsystem_type) :: velo_grid !> coordinate system of the velocity grid

    real(dp), dimension(:),pointer :: x0 => null() !original x0 grid 
    real(dp), dimension(:),pointer :: y0 => null() !original y0 grid
    real(dp), dimension(:),pointer :: x1 => null() !original x1 grid
    real(dp), dimension(:),pointer :: y1 => null() !original y1 grid

    integer :: global_bc = 0     ! 0 for periodic, 1 for outflow, 2 for no-penetration

  end type glide_general

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_options

    !> Holds user options controlling the methods used in the ice-model integration.

    !-----------------------------------------------------------------------
    ! standard options
    !-----------------------------------------------------------------------

    integer :: whichdycore = 2

    ! Choice of two dycores:
    !> \begin{description} 
    !> \item[0] Glide dycore (SIA, serial (SLAP) only)
    !> \item[1] Glam dycore (no longer supported!)
    !> \item[2] Glissade dycore (1st-order, FEM, serial (SLAP) or parallel (F90 native PCG solver) )
    !> \item[3] FELIX-Albany dycore (1st-order, FEM, using Trilino/Albany, mesh information from Glissade)
    !> \item[4] BISICLES dycore (L1L2, FVM, parallel using Chombo AMR)
    !> \end{description}

    integer :: whichevol = 0

    !TODO: Break into two options: one set (0, 1, 2) for Glide and one set (3, 4, 5) for Glissade
    !> Thickness evolution method:
    !> \begin{description}
    !> \item[0] Pseudo-diffusion 
    !> \item[1] Alternating direction implicit (ADI)
    !> \item[2] Diffusion (also calculates velocities) 
    !> \item[3] Incremental remapping
    !> \item[4] 1st-order upwind scheme
    !> \item[5] Temperature evolves but thickness does not
    !> \end{description}

    integer :: whichtemp = 1

    !> Method of ice temperature calculation:
    !> \begin{description} 
    !> \item[0] Set column to surface air temperature
    !> \item[1] Prognostic temperature solution 
    !> \item[2] Do NOTHING - hold temperatures steady at initial value  
    !> \item[3] Prognostic enthalpy solution
    !> \end{description}

    integer :: temp_init = 1

    ! Temperature initialization:
    !> \begin{description} 
    !> \item[0] Initialize temperature to 0 C
    !> \item[1] Initialize temperature to surface air temperature
    !> \item[2] Initialize temperature with a linear profile in each column
    !> \item[3] Initialize temperature with an advective-diffusive balance in each column
    !> \item[4] Initialize temperature from external file
    !> \end{description}

    !> Method for calculating flow factor $A$:

    integer :: whichflwa = 2

    !> \begin{description} 
    !> \item[0] Set equal to $1\times 10^{-16} \mathrm{Pa}^{-n} \mathrm{yr}^{-1}$
    !> \item[1] \emph{Paterson and Budd} relationship, 
    !> with temperature set to $-5^{\circ}\mathrm{C}$ 
    !> \item[2] \emph{Paterson and Budd} relationship
    !> \item[3] Read flwa/flwastag from file
    !> \end{description}

    integer :: whichbtrc = 0

    !> Basal slip coefficient (used for glide and for glissade local SIA; not glissade HO):
    !> \begin{description}
    !> \item[0] Set equal to zero everywhere
    !> \item[1] Set to (non--zero) constant
    !> \item[2] Set to (non--zero) constant where basal water is present, otherwise to zero
    !> \item[3] Set to (non--zero) constant where bed temperature is at pressure melting point, otherwise to zero
    !> \item[4] linear function of basal melt rate
    !> \item[5] \texttt{tanh} function of basal water depth 
    !> \end{description}

    integer :: whichbwat = 0

    !> Basal water depth: 
    !> \begin{description} 
    !> \item[0] Set to zero everywhere 
    !> \item[1] Compute from local basal water balance 
    !> \item[2] Compute the basal water flux, then find depth via calculation
    !> \item[3] Set to constant (10 m) everywhere, to force T = Tpmp.
    !> \end{description}

    integer :: whichbmlt_float = 0

    !> basal melt rate for floating ice:
    !> \begin{description}
    !> \item[0] Basal melt rate = 0 for floating ice
    !> \item[1] Depth-dependent basal melt rate for floating ice as specified for MISMIP+
    !> \item[2] Basal melt rate = constant for floating ice (with option to selectively mask out melting)
    !> \item[3] Depth-dependent basal melt rate for floating ice
    !> \item[4] External basal melt rate field (from input file or coupler)
    !> \item[5] Basal melt rate for floating ice from MISOMIP ocean forcing with plume model
    !> \end{description}

    logical :: enable_bmlt_anomaly = .false.
    !> if true, then apply a prescribed anomaly to bmlt_float

    integer :: basal_mbal = 1

    !> basal mass balance:
    !> \begin{description}
    !> \item[0] Basal mass balance not included in continuity equation
    !> \item[1] Basal mass balance included in continuity equation
    !> \end{description}

    integer :: smb_input = 0

    !> units for SMB input:
    !> \begin{description}
    !> \item[0] SMB input in units of m/yr ice (same as acab)
    !> \item[1] SMB input in units of mm/yr water equivalent
    !> \end{description}
    
    logical :: enable_acab_anomaly = .false.
    !> if true, then apply a prescribed anomaly to acab

    integer :: overwrite_acab = 0
    !> overwrite acab (m/yr ice) in selected regions:
    !> \begin{description}
    !> \item[0] Do not overwrite acab anywhere
    !> \item[1] Overwrite acab where input acab = 0
    !> \item[2] Overwrite acab where input thickness <= threshold value
    !> \end{description}

    integer :: gthf = 0

    !> geothermal heat flux:
    !> \begin{description}
    !> \item[0] prescribed uniform geothermal heat flux
    !> \item[1] read 2D geothermal flux field from input file (if present)
    !> \item[2] calculate geothermal flux using 3d diffusion
    !> \end{description}

    ! This replaces model%isos%do_isos
    integer :: isostasy = 0

    !> isostasy:
    !> \begin{description}
    !> \item[0] no isostatic adjustment
    !> \item[1] compute isostatic adjustment using lithosphere/asthenosphere model
    !> \end{description}

    integer :: whichcalving = 1

    !> Calving: 
    !> \begin{description} 
    !> \item[0] No calving 
    !> \item[1] Set thickness to zero if floating
    !> \item[2] Lose a fraction of floating ice at marine margin
    !> \item[3] Set thickness to zero if relaxed bedrock is more than a
    !>          certain water depth (variable "marine_limit" in glide_types)  
    !> \item[4] Set thickness to zero if present bedrock topography lies below
    !>          a certain water depth (variable "marine_limit" in glide_types)  
    !> \item[5] Set thickness to zero based on grid location (field 'calving_mask')
    !> \item[6] Set thickness to zero if ice at marine margin is thinner than
    !>          a certain value (variable 'calving_minthck' in glide_types)
    !> \item[7] Set thickness to zero based on stress (eigencalving) criterion
    !> \item[8] Calve ice that is sufficiently damaged
    !> \item[9] Huybrechts grounding line scheme for Greenland initialization
    !> \end{description}

    integer :: calving_init = 0
    !> \begin{description}
    !> \item[0] Do not calve at initialization
    !> \item[1] Calve at initialization
    !> \end{description}

    integer :: calving_domain = 1
    !> \begin{description}
    !> \item[0] Calve only at ocean edge
    !> \item[1] Calve wherever the calving criterion is met
    !> \end{description}

    logical  :: remove_icebergs = .true.
    !> if true, then identify and remove icebergs after calving
    !> These are connected regions with zero basal traction and no connection to grounded ice.
    !>       Safer to make it true, but not necessary for all applications

    logical  :: limit_marine_cliffs = .false.
    !> if true, then thin marine-based cliffs based on a thickness threshold

    logical  :: cull_calving_front = .false.
    !> if true, then cull calving_front cells at initialization
    !> This can make the run more stable by removing long, thin peninsulas

    integer :: dm_dt_diag = 0
    !> \begin{description}
    !> \item[0] Write dmass/dt diagnostic in units of kg/s
    !> \item[1] Write dmass/dt diagnostic in units of Gt/yr
    !> \end{description}

    integer :: diag_minthck = 1
    !> \begin{description}
    !> \item[0] Include cells with H > 0 in global diagnostics
    !> \item[1] Include cells with H > thklim in global diagnostics
    !> \end{description}

    integer :: whichwvel = 0

    !> Vertical velocities: 
    !> \begin{description}
    !> \item[0] Usual vertical integration 
    !> \item[1] Vertical integration constrained to obey upper kinematic BC
    !> \end{description}

    integer :: which_sigma = 0

    !> \begin{description}
    !> \item[0] compute standard Glimmer sigma coordinates
    !> \item[1] sigma coordinates are given in external file
    !> \item[2] sigma coordinates are given in configuration file
    !> \item[3] evenly spaced levels
    !> \item[4] compute Pattyn sigma coordinates
    !> \end{description}

    !TODO - Make is_restart a logical variable?
    integer :: is_restart = 0
    !> if the run is a restart of a previous run
    !> \begin{description}
    !> \item[0] normal start-up (take init fields from .nc input file OR if absent, use default options)
    !> \item[1] restart model from previous run (do not calc. temp, rate factor, or vel)
    !> \end{description}

    integer :: restart_extend_velo = 0
    !> if velocity fields should be written on the extended staggered mesh
    !> \begin{description}
    !> \item[0] write uvel and vvel to restart file on standard staggered mesh
    !> \item[1] write uvel_extend and vvel_extend to restart file on extended staggered mesh
    !>          (required if restart velocities are nonzero on global boundaries)
    !> \end{description}

    ! This is a Glimmer serial option
    ! The parallel code enforces periodic EW and NS boundary conditions by default
    logical :: periodic_ew = .false.

    !> \begin{description}
    !> \item[0] no periodic EW boundary conditions
    !> \item[1] periodic EW boundary conditions
    !> \end{description}

    !-----------------------------------------------------------------------
    ! Higher-order options associated with the glissade dycore
    !-----------------------------------------------------------------------

    integer :: which_ho_efvs = 2

    !> Flag that indicates how effective viscosity is computed
    !> \begin{description}
    !> \item[0] constant value
    !> \item[1] multiple of flow factor
    !> \item[2] compute from effective strain rate

    integer :: which_ho_disp = 1

    !> Flag that indicates method for computing the dissipation during the temperature calc.
    !> \begin{description}
    !> \item[-1] for no dissipation
    !> \item[0] for 0-order SIA approx
    !> \item[1] for first-order dissipation (Blatter-Pattyn)
    !>      
    !> \end{description}

    integer :: which_ho_thermal_timestep = 0
    !> \begin{description}
    !> \item[0] vertical thermal solve before transport solve
    !> \item[1] vertical thermal solve after transport solve
    !> \item[2] vertical thermal solve split; both before and after transport solve
    !> \end{description}

    integer :: which_ho_babc = 4

    !> Flag that describes basal boundary condition for HO dyn core: 
    !> \begin{description}
    !> \item[0] spatially uniform value (low value of 10 Pa/yr by default)
    !> \item[1] large value for frozen bed, lower value for bed at pressure melting point
    !> \item[2] treat beta value as a till yield stress (in Pa) using Picard iteration 
    !> \item[3] pseudo-plastic basal sliding law; can model linear, power-law or plastic behavior
    !> \item[4] very large value for beta to enforce no slip everywhere 
    !> \item[5] beta field passed in from .nc input file as part of standard i/o
    !> \item[6] no slip everywhere (using Dirichlet BC rather than large beta)
    !> \item[7] treat beta value as till yield stress (in Pa) using Newton-type iteration (in development)
    !> \item[8] beta field as prescribed for ISMIP-HOM test C (serial only)
    !> \item[9] power law
    !> \item[10] Coulomb friction law using effective pressure, with flwa from lowest ice layer
    !> \item[11] Coulomb friction law using effective pressure, with constant basal flwa
    !> \item[12] basal stress is the minimum of Coulomb and power-law values, as in Tsai et al. (2015)
    !> \item[13] power law using effective pressure
    !> \item[14] simple hard-coded pattern (useful for debugging)
    !> \end{description}

    integer :: which_ho_inversion = 0
    !> Flag for basal traction inversion options
    !> Note: Inversion is currently supported for which_ho_babc = 11 only
    !> \begin{description}
    !> \item[0] no inversion
    !> \item[1] invert for basal traction parameters and subshelf melting
    !> \item[2] apply parameters and subshelf melting from previous inversion
    !> \end{description}

    integer :: which_ho_bwat = 0
    !> Basal water depth:
    !> \begin{description}
    !> \item[0] Set to zero everywhere
    !> \item[1] Set to constant everywhere, to force T = Tpmp.
    !> \item[2] Local basal till model with constant drainage
    !> \end{description}

    integer :: which_ho_effecpress = 0
    !> Flag that describes effective pressure calculation for HO dyn core: 
    !> \begin{description}
    !> \item[0] N = overburden pressure, rhoi*grav*thck
    !> \item[1] N is reduced where the bed is at or near the pressure melting point
    !> \item[2] N is reduced where there is melting at the bed
    !> \item[3] N is reduced due to connection of subglacial water to the ocean
    !> \item[4] N is reduced where basal water is present
    !> \end{description}

    integer :: which_ho_nonlinear = 0
    !> Flag that indicates method for solving the nonlinear iteration when solving 
    !> the first-order momentum balance
    !> \begin{description}
    !> \item[0] use the standard Picard iteration
    !> \item[1] use Jacobian Free Newton Krylov (JFNK) method
    !> \end{description}

    integer :: which_ho_resid = 3
    !> Flag that indicates method for computing residual in PP dyn core: 
    !> \begin{description}
    !> \item[0] maxval 
    !> \item[1] maxval ignoring basal velocity 
    !> \item[2] mean value
    !> \item[3] L2 norm of system residual, Ax-b=resid
    !> \item[4] L2 norm of system residual relative to rhs, |Ax-b|/|b|
    !> \end{description}

    integer :: which_ho_sparse = 0
    !> Flag that indicates method for solving the sparse linear system
    !> that arises from the higher-order solver
    !> \begin{description}
    !> \item[-1] SLAP (serial): Preconditioned conjugate gradient, incomplete Cholesky preconditioner
    !> \item[0]  SLAP (serial): Biconjugate gradient, incomplete LU preconditioner
    !> \item[1]  SLAP (serial): GMRES, incomplete LU preconditioner
    !> \item[2] Native PCG, parallel-enabled, standard solver
    !> \item[3] Native PCG, parallel-enabled, Chronopoulos-Gear solver
    !> \item[4] standalone interface to Trilinos
    !> \end{description}

    ! parameters to store external dycore options/information -- Doug Ranken 04/20/12
    integer*4 :: external_dycore_type = 0
    integer*4 :: external_dycore_model_index = -1  
    !> Flag to select an external dynamic core.
    !> \begin{description}
    !> \item[0] Do not use an external dynamic core
    !> \item[1] Use the BISICLES external dynamic core
    !> \item[2] Use the ALBANY_FELIX external dynamic core
    !> \end{description}

    character(fname_length) :: dycore_input_file=''
    !FD Name of a file containing external dycore settings.

    integer :: which_ho_approx = 2 
    !> Flag that indicates which Stokes approximation to use with the glissade dycore.
    !> Not valid for other dycores 
    !> Compute Blatter-Pattyn HO momentum balance by default.
    !> Note: There are two SIA options:
    !>       Option -1 uses module glissade_velo_sia to compute local SIA velocities, similar to Glide
    !>       Option 0 uses module glissade_velo_higher to compute SIA velocities via an iterative solve
    !> \begin{description}
    !> \item[-1] Shallow-ice approximation, Glide-type calculation; uses glissade_velo_sia
    !> \item[0]  Shallow-ice approximation, vertical-shear stresses only; uses glissade_velo_higher
    !> \item[1]  Shallow-shelf approximation, horizontal-plane stresses only; uses glissade_velo_higher
    !> \item[2]  Blatter-Pattyn approximation with both vertical-shear and horizontal-plane stresses; uses glissade_velo_higher
    !> \item[3]  Vertically integrated 'L1L2' approximation with vertical-shear and horizontal-plane stresses; uses glissade_velo_higher
    !> \item[4]  Depth-integrated viscosity approximation based on Goldberg (2011); uses glissade_velo_higher 
    !> \end{description}

    integer :: which_ho_precond = 2    
    !> Flag that indicates which Stokes preconditioner to use in the glissade dycore.
    !> Not valid for other dycores 
    !> \begin{description}
    !> \item[0] No preconditioner
    !> \item[1] Diagonal preconditioner
    !> \item[2] Physics-based shallow-ice preconditioner
    !> \end{description}

    integer :: which_ho_gradient = 0    
    !> Flag that indicates which gradient operator to use in the glissade dycore.
    !> Not valid for other dycores
    !> NOTE: Upstream may be better for ice evolution because it damps checkerboard noise.
    !> \begin{description}
    !> \item[0] Centered gradient
    !> \item[1] First-order accurate upstream gradient
    !> \item[2] Second-order accurate upstream gradient

    !WHL - Changed default from 1 to 3.
    !      Option 3 is appropriate for ice sheets with both land and marine boundaries,
    !       when a lateral spreading force is computed for marine ice cliffs.
    !      This lateral force ensures spreading where the gradient is zero.

    integer :: which_ho_gradient_margin = 1
    !> Flag that indicates how to compute the gradient at the ice margin in the glissade dycore.
    !> Note: Gradients are always computed at edges with ice on both sides.
    !>       The methods differ in whether gradients are computed when an ice-covered cell
    !>        lies above an ice-free cell (land or ocean).
    !> Not valid for other dycores
    !> \begin{description}
    !> \item[0] Compute edge gradient when either cell is ice-covered
    !> \item[1] Compute edge gradient for ice-covered cell above ice-free land (not ocean)
    !> \item[2] Compute edge gradient only when both cells have ice

    !TODO: Change the default to 2nd order vertical remapping
    ! WHL: Keeping 1st order vertical remapping for now, pending more testing
    integer :: which_ho_vertical_remap = 0
    !> Flag that indicates the order of accuracy for vertical remapping
    !> \begin{description}
    !> \item[0] first-order accurate in the vertical direction
    !> \item[1] second-order accurate in the vertical direction

    integer :: which_ho_assemble_taud = 1

    !> Flag that describes how driving-stress terms are assembled in the glissade finite-element calculation
    !> \begin{description}
    !> \item[0] standard finite-element calculation (which effectively smooths the driving stress)
    !> \item[1] apply local value of driving stress at each vertex

    integer :: which_ho_assemble_beta = 1

    !> Flag that describes how beta terms are assembled in the glissade finite-element calculation
    !> \begin{description}
    !> \item[0] standard finite-element calculation (which effectively smooths beta at discontinuities)
    !> \item[1] apply local value of beta at each vertex

    integer :: which_ho_assemble_bfric = 1

    !> Flag that describes how the basal friction heat flux is computed in the glissade finite-element calculation
    !> \begin{description}
    !> \item[0] standard finite-element calculation summing over quadrature points
    !> \item[1] apply local value of beta*(u^2 + v^2) at each vertex

    integer :: which_ho_calving_front = 0
    !> Flag that indicates whether to use a subgrid calving front parameterization
    !> \begin{description}
    !> \item[0] no subgrid calving front parameterization
    !> \item[1] subgrid calving front parameterization

    integer :: which_ho_ground = 0
    !> Flag that indicates how to compute the grounded fraction of each gridcell in the glissade dycore.
    !> Not valid for other dycores
    !> \begin{description}
    !> \item[0] fground = 0 in floating cells (based on flotation condition), else fground = 1 
    !> \item[1] 0 <= fground <= 1, based on a grounding line parameterization
    !> \item[2] similar to [1], but f_ground is summed over quadrants rather than staggered cells

    integer :: which_ho_ground_bmlt = 0
    !> Flag that indicates how to compute bmlt_float in partly grounded cells
    !> \begin{description}
    !> \item[0] Apply bmlt_float in all floating cells, based on floating_mask
    !> \item[1] Weigh bmlt by floating fraction of cell if 0 < f_ground < 1
    !> \item[2] Set bmlt_float = 0 in partly grounded cells (f_ground > 0)

    integer :: which_ho_flotation_function = 2
    !> Flag that indicates how to compute the flotation function at and near vertices in the glissade dycore
    !> Not valid for other dycores
    !> \begin{description}
    !> \item[0] f_flotation = (-rhow*b/rhoi*H) = f_pattyn; <=1 for grounded, > 1 for floating
    !> \item[1] f_flotation = (rhoi*H)/(-rhow*b) = 1/f_pattyn; >=1 for grounded, < 1 for floating
    !> \item[2] f_flotation = -rhow*b - rhoi*H = ocean cavity thickness; <=0 for grounded, > 0 for floating 

    logical :: block_inception = .false.
    !> Flag that indicates whether ice inception away from the main ice sheet is blocked

    logical :: remove_ice_caps = .false.
    !> Flag that indicates whether ice caps are removed and added to the calving flux

    integer :: which_ho_ice_age = 1
    !> Flag that indicates whether to compute a 3d ice age tracer
    !> \item[0] ice age computation off
    !> \item[1] ice age computation on

    integer :: glissade_maxiter = 100    
    !> maximum number of nonlinear iterations to be used by the Glissade velocity solver

    ! The remaining options are not currently supported

    !integer :: which_bproc = 0
    !Options for the basal processes code
    !> \begin{description}
    !> \item[0] Disabled
    !> \item[1] Full calculation, with at least 3 nodes to represent the till layer
    !> \item[2] Fast calculation, using Tulaczyk empirical parametrization
    !> \end{description}


  end type glide_options

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_geometry

    !> Holds fields and other information relating to the
    !> geometry of the ice sheet and bedrock.

    real(dp),dimension(:,:),pointer :: thck => null()
    !> The thickness of the ice, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: usrf => null()
    !> The elevation of the upper ice surface, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: lsrf => null() 
    !> The elevation of the lower ice surface, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: topg => null() 
    !> The elevation of the topography, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: thck_obs => null()
    !> Observed ice thickness, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: usrf_obs => null()
    !> Observed upper surface elevation, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: topg_obs => null()
    !> Observed basal topography, divided by \texttt{thk0}.

    real(dp),dimension(:,:),pointer :: f_flotation => null() 
    !> flotation function, > 0 for floating ice and <= 0 for grounded ice

    real(dp),dimension(:,:),pointer :: f_ground => null() 
    !> The fractional area at each vertex which is grounded

    real(dp),dimension(:,:),pointer :: f_ground_cell => null() 
    !> The fractional area in each cell which is grounded

    real(dp),dimension(:,:,:),pointer :: ice_age => null()
    !> The age of a given ice layer, divided by \texttt{tim0}.
    !> Used to be called 'age', but changed to 'ice_age' for easier grepping

    real(dp),dimension(:,:),pointer :: thck_old => null()        !> old ice thickness, divided by \texttt{thk0}
    real(dp),dimension(:,:),pointer :: dthck_dt => null()        !> ice thickness tendency (m/s)
    real(dp),dimension(:,:),pointer :: dthck_dt_tavg => null()   !> ice thickness tendency (m/s, time average)

    real(dp),dimension(:,:),pointer :: cell_area => null()
    !> The cell area of the grid, divided by \texttt{len0*len0}.

    integer :: ntracers
    !> number of tracers to be transported

    real(dp),dimension(:,:,:,:), pointer :: tracers => null()
    !> all 3D tracers packed into one array

    real(dp),dimension(:,:,:), pointer :: tracers_usrf => null()
    !> value at upper surface for all tracers

    real(dp),dimension(:,:,:), pointer :: tracers_lsrf => null()
    !> value at lower surface for all tracers

    integer, dimension(:,:),pointer :: thkmask => null()
    !> see glide_mask.f90 for possible values

    integer, dimension(:,:),pointer :: stagmask => null()
    !> see glide_mask.f90 for possible values

    ! mass fluxes at upper, lower and lateral boundaries
    ! TODO: Move to a flux derived type?
    ! Note: sfc_mbal_flux and basal_mbal_flux are not strictly needed, since they are equal to acab_applied and bmlt_applied
    !       multipled by a constant. For some applications, however, it may be useful to output the mass balance in SI units.
    real(dp),dimension(:,:),  pointer :: sfc_mbal_flux =>null()        !> surface mass balance (kg m^-2 s^-1), diagnosed from acab
    real(dp),dimension(:,:),  pointer :: sfc_mbal_flux_tavg =>null()   !> surface mass balance (kg m^-2 s^-1, time average)
    real(dp),dimension(:,:),  pointer :: basal_mbal_flux =>null()      !> basal mass balance (kg m^-2 s^-1), diagnosed from bmlt
    real(dp),dimension(:,:),  pointer :: basal_mbal_flux_tavg =>null() !> basal mass balance (kg m^-2 s^-1, time average)
    real(dp),dimension(:,:),  pointer :: calving_flux =>null()         !> calving flux (kg m^-2 s^-1), diagnosed from calving_thck
    real(dp),dimension(:,:),  pointer :: calving_flux_tavg =>null()    !> calving flux (kg m^-2 s^-1, time average)
    real(dp),dimension(:,:),  pointer :: gl_flux_east =>null()         !> mass flux eastward at grounding line, edge-based (kg m^-1 s^-1)
    real(dp),dimension(:,:),  pointer :: gl_flux_north =>null()        !> mass flux northward at grounding line, edge_based (kg m^-1 s^-1)
    real(dp),dimension(:,:),  pointer :: gl_flux =>null()              !> mass flux at grounding line, cell-based (kg m^-1 s^-1)
    real(dp),dimension(:,:),  pointer :: gl_flux_tavg =>null()         !> mass flux at grounding line, cell-based (kg m^-1 s^-1, time average)

    !* (DFM ----------------- The following fields were added for BISICLES interface --------------)
    !*SFP: These fields need to be passed to POP for ice ocean coupling
    ! WHL: When Dan added the masks, he made them real-valued. They are now integers, which might break the POP coupling.
    real(dp),dimension(:,:),pointer :: lower_cell_loc => null()  !> z-location of the center of the lowest ice cell center
    real(dp),dimension(:,:),pointer :: lower_cell_temp => null() !> temperature in the cell located at lower_cell_loc
    integer, dimension(:,:),pointer :: ice_mask => null()        !> = 1 where ice is present, else = 0
    integer, dimension(:,:),pointer :: ice_mask_stag => null()   !> = 1 where ice is present on staggered grid, else = 0
    integer, dimension(:,:),pointer :: floating_mask => null()   !> = 1 where ice is present and floating, else = 0
    integer, dimension(:,:),pointer :: grounded_mask => null()   !> = 1 where ice is present and grounded, else = 0
    integer, dimension(:,:),pointer :: ice_sheet_mask => null()  !> = 1 for ice sheet cells, = 0 for ice cap cells
    integer, dimension(:,:),pointer :: ice_cap_mask => null()    !> = 1 for ice cap cells disconnected from the main ice sheet

    integer, dimension(:,:),pointer :: thck_index => null()
    ! Set to nonzero integer for ice-covered cells (thck > 0), cells adjacent to ice-covered cells,
    !  and cells with acab > 0.  The non-zero points are numbered in sequence from the bottom left 
    !  to the top right, going along the rows.

    integer :: totpts = 0       ! total number of points with nonzero thck_index
    logical :: empty = .true.   ! true if totpts = 0

    ! global scalars

    real(dp) :: iarea                  ! total ice area (m^2)
    real(dp) :: iareag                 ! total grounded ice area (m^2)
    real(dp) :: iareaf                 ! total floating ice area (m^2)
    real(dp) :: ivol                   ! total ice volume (m^3)
    real(dp) :: imass                  ! total ice mass (kg)
    real(dp) :: imass_above_flotation  ! total ice mass above flotation (kg)
    real(dp) :: total_smb_flux         ! total surface mass balance flux (kg/s)
    real(dp) :: total_bmb_flux         ! total basal mass balance flux (kg/s)
    real(dp) :: total_calving_flux     ! total calving mass flux (kg/s)
    real(dp) :: total_gl_flux          ! total grounding line mass flux (kg/s)
    real(dp) :: total_smb_flux_tavg    ! total surface mass balance flux (kg/s), time average
    real(dp) :: total_bmb_flux_tavg    ! total basal mass balance flux (kg/s), time average
    real(dp) :: total_calving_flux_tavg! total calving mass flux (kg/s), time average
    real(dp) :: total_gl_flux_tavg     ! total grounding line mass flux (kg/s), time average

  end type glide_geometry

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_geomderv

    !> Holds the horizontal and temporal derivatives of the thickness and
    !> upper surface elevation, as well as the thickness on the staggered grid.

    !*tb* Added a bunch of stuff here to clean up the higher order code that
    !I've been writing.  Might be worth it to add a mechanism to conditionally
    !allocate these depending on whether they are needed by the SIA core or by
    !the higher-order extensions

    !First derivatives on a staggered grid
    real(dp),dimension(:,:),pointer :: dthckdew => null() !> E-W derivative of thickness.
    real(dp),dimension(:,:),pointer :: dusrfdew => null() !> E-W derivative of upper surface elevation.
    real(dp),dimension(:,:),pointer :: dthckdns => null() !> N-S derivative of thickness.
    real(dp),dimension(:,:),pointer :: dusrfdns => null() !> N-S derivative of upper surface elevation.
    real(dp),dimension(:,:),pointer :: dlsrfdew => null() !*tb* added
    real(dp),dimension(:,:),pointer :: dlsrfdns => null() !*tb* added

    !Time derivatives
    real(dp),dimension(:,:),pointer :: dthckdtm => null() !> Temporal derivative of thickness.
    real(dp),dimension(:,:),pointer :: dusrfdtm => null() !> Temporal derivative of upper surface elevation.

    !TODO - Move staggered variables from glide_geomderv type to glide_geometry?

    !Staggered grid versions of geometry variables
    real(dp),dimension(:,:),pointer :: stagthck => null() !> Thickness averaged onto the staggered grid.
    real(dp),dimension(:,:),pointer :: stagusrf => null() !> Upper surface averaged onto the staggered grid
    real(dp),dimension(:,:),pointer :: staglsrf => null() !> Lower surface averaged onto the staggered grid
    real(dp),dimension(:,:),pointer :: stagtopg => null() !> Bedrock topography averaged onto the staggered grid

  end type glide_geomderv

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_tensor
    real(dp), dimension(:,:,:), pointer :: scalar => null()
    real(dp), dimension(:,:,:), pointer :: xz => null()
    real(dp), dimension(:,:,:), pointer :: yz => null()
    real(dp), dimension(:,:,:), pointer :: xx => null()
    real(dp), dimension(:,:,:), pointer :: yy => null()
    real(dp), dimension(:,:,:), pointer :: xy => null()
  end type glide_tensor
  
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_velocity

    !> Holds the velocity fields in 2D and 3D. At least some of these fields
    real(dp),dimension(:,:,:),pointer :: uvel  => null()   !> 3D $x$-velocity.
    real(dp),dimension(:,:,:),pointer :: vvel  => null()   !> 3D $y$-velocity.
    real(dp),dimension(:,:,:),pointer :: velnorm => null() ! horizontal ice speed
    real(dp),dimension(:,:,:),pointer :: wvel  => null()   !> 3D $z$-velocity.
    real(dp),dimension(:,:,:),pointer :: wgrd  => null()   !> 3D grid vertical velocity.
    real(dp),dimension(:,:)  ,pointer :: uflx  => null()   !> 
    real(dp),dimension(:,:)  ,pointer :: vflx  => null()   !> 
    real(dp),dimension(:,:)  ,pointer :: diffu => null()   !> 
    real(dp),dimension(:,:)  ,pointer :: diffu_x => null() !*sfp* moved from velocity_hom deriv type
    real(dp),dimension(:,:)  ,pointer :: diffu_y => null() 
    real(dp),dimension(:,:)  ,pointer :: total_diffu => null() !> total diffusivity

    ! Note: DIVA solves for uvel_2d and vvel_2d; these are typically (but not necessarily) the vertical average
    real(dp),dimension(:,:)  ,pointer :: uvel_2d  => null()   !> 2D $x$-velocity; typically the vertical average
    real(dp),dimension(:,:)  ,pointer :: vvel_2d  => null()   !> 2D $y$-velocity; typically the vertical average
    real(dp),dimension(:,:)  ,pointer :: ubas  => null()      !> basal $x$-velocity
    real(dp),dimension(:,:)  ,pointer :: vbas  => null()      !> basal $y$-velocity
    real(dp),dimension(:,:)  ,pointer :: uvel_mean  => null() !> vertical mean $x$-velocity
    real(dp),dimension(:,:)  ,pointer :: vvel_mean  => null() !> vertical mean $y$-velocity

    ! Note: uvel_extend and vvel_extend can be used for input and output of uvel, vvel on a staggered grid 
    !       that is the same size as the unstaggered grid. This is required for exact restart if velocities
    !       are nonzero along the north and east boundaries of the global domain.
    real(dp),dimension(:,:,:),pointer :: uvel_extend => null()  !> 3D $x$-velocity on extended staggered grid
    real(dp),dimension(:,:,:),pointer :: vvel_extend => null()  !> 3D $y$-velocity on extended staggered grid
    real(dp),dimension(:,:)  ,pointer :: uvel_2d_extend => null()  !> 2D $x$-velocity on extended staggered grid
    real(dp),dimension(:,:)  ,pointer :: vvel_2d_extend => null()  !> 2D $y$-velocity on extended staggered grid

    real(dp),dimension(:,:)  ,pointer :: bed_softness => null() !> bed softness parameter
    real(dp),dimension(:,:)  ,pointer :: btrc  => null()        !>  basal traction (scaler field)
    real(dp),dimension(:,:,:),pointer :: btraction => null()    !> x(1,:,:) and y(2,:,:) "consistent" basal traction fields 
    real(dp),dimension(:,:)  ,pointer :: beta  => null()        !> basal shear coefficient on velo grid (Pa yr/m)
    real(dp),dimension(:,:)  ,pointer :: beta_internal => null()!> beta weighted by f_ground or otherwise adjusted (glissade only)
    real(dp),dimension(:,:)  ,pointer :: unstagbeta  => null()  !> basal shear coefficient on ice grid (Pa yr/m)
    real(dp),dimension(:,:)  ,pointer :: tau_x => null()        !> SIA basal shear stress, x-dir
    real(dp),dimension(:,:)  ,pointer :: tau_y => null()        !> SIA basal shear stress, y-dir

    !> mask that specifies where the velocity being read in should be held constant as a dirichlet condition
    integer, dimension(:,:), pointer  :: kinbcmask => null()    

    !> masks that specify where the outflow velocities on the global boundary should be set to zero
    integer, dimension(:,:), pointer  :: umask_no_penetration => null()
    integer, dimension(:,:), pointer  :: vmask_no_penetration => null()

    ! for viewing the spatial pattern of residuals
    real(dp),dimension(:,:,:),pointer :: resid_u => null()     ! u component of residual Ax - b where x is the velocity
    real(dp),dimension(:,:,:),pointer :: resid_v => null()     ! v component of residual Ax - b where x is the velocity

    ! for viewing the driving stress on the RHS
    real(dp),dimension(:,:,:),pointer :: rhs_u => null()       ! u component of b in Ax = b
    real(dp),dimension(:,:,:),pointer :: rhs_v => null()       ! v component of b in Ax = b

    ! strain rate tensor (HO only) 
    type(glide_tensor) :: strain_rate                          ! strain rate tensor, diagnosed from stress tensor and efvs (s^-1)

    real(dp), dimension(:,:), pointer :: divu => null()        ! horizontal divergence rate, eps_xx + eps_yy (s^-1)

  end type glide_velocity

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_stress_t      

    type(glide_tensor) :: tau         ! HO only
    real(dp),dimension(:,:,:),pointer :: efvs => null()    !> effective viscosity
    real(dp),dimension(:,:),  pointer :: efvs_vertavg => null()  !> vertical mean effective viscosity
    real(dp),dimension(:,:),  pointer :: btractx => null() !> basal traction (Pa), x comp
    real(dp),dimension(:,:),  pointer :: btracty => null() !> basal traction (Pa), y comp
    real(dp),dimension(:,:),  pointer :: btract => null()  !> basal traction (Pa), magnitude = sqrt(btractx^2 + btracty^2)
    !WHL - The extended versions are needed for exact restart if using DIVA solver for a problem with nonzero traction at global boundaries
    real(dp),dimension(:,:),  pointer :: btractx_extend => null() !> basal traction (Pa), x comp, on extended staggered grid
    real(dp),dimension(:,:),  pointer :: btracty_extend => null() !> basal traction (Pa), y comp, on extended staggered grid
    real(dp),dimension(:,:),  pointer :: taudx => null()   !> driving stress (Pa), x comp
    real(dp),dimension(:,:),  pointer :: taudy => null()   !> driving stress (Pa), y comp

  end type glide_stress_t      

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!TODO - Make eus a config file parameter.
!TODO - Rename acab in glide_climate type to avoid confusion over units?
!       Here, acab has units of m/y ice, whereas in Glint, acab has units of m/y water equiv.

  ! Note on acab_tavg: This is the average value of acab over an output interval.
  !                    If 'average = 1' in the acab entry of glide_vars.def, then acab_tavg is automatically
  !                     accumulated and averaged during runtime, without any additional code needed.
  !                    Other variables with a '_tavg' suffix are handled similarly.
  !
  ! Note on acab_corrected: Optionally, acab can be supplemented with a flux correction or an anomaly.
  !                         The background field, acab, does not include the corrections.
  !                         Write acab_corrected to the output file to see the modified SMB field.

  type glide_climate
     !> Holds fields used to drive the model
     real(dp),dimension(:,:),pointer :: acab            => null() !> Surface mass balance (m/yr ice)
     real(dp),dimension(:,:),pointer :: acab_tavg       => null() !> Surface mass balance (time average).
     real(dp),dimension(:,:),pointer :: acab_anomaly    => null() !> Surface mass balance anomaly (m/yr ice)
     real(dp),dimension(:,:),pointer :: acab_corrected  => null() !> Surface mass balance with flux or anomaly corrections (m/yr ice)
     real(dp),dimension(:,:),pointer :: acab_applied    => null() !> Surface mass balance applied to ice (m/yr ice)
                                                                  !>    = 0 for ice-free cells with acab < 0
     real(dp),dimension(:,:),pointer :: acab_applied_tavg => null() !> Surface mass balance applied to ice (m/yr ice, time average)
     real(dp),dimension(:,:),pointer :: smb             => null() !> Surface mass balance (mm/yr water equivalent)
                                                                  !> Note: acab (m/y ice) is used internally by dycore, 
                                                                  !>       but can use smb (mm/yr w.e.) for I/O
     real(dp),dimension(:,:),pointer :: artm            => null() !> Annual mean air temperature (degC)
     integer, dimension(:,:),pointer :: no_advance_mask => null() !> mask for cells where advance is not allowed 
                                                                  !> (any ice reaching these locations is eliminated)
     integer, dimension(:,:),pointer :: overwrite_acab_mask => null() !> mask for cells where acab is overwritten

     real(dp) :: eus = 0.d0                         !> eustatic sea level
     real(dp) :: acab_factor = 1.0d0                !> adjustment factor for external acab field (unitless)
     real(dp) :: acab_anomaly_timescale = 0.0d0     !> number of years over which the acab anomaly is phased in linearly
                                                    !> If set to zero, then the anomaly is applied immediately.
                                                    !> The initMIP value is 40 yr.
     real(dp) :: overwrite_acab_value = 0.0d0       !> acab value to apply in grid cells where overwrite_acab_mask = 1
     real(dp) :: overwrite_acab_minthck = 0.0d0     !> overwrite acab where thck <= overwrite_acab_minthck

  end type glide_climate

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_calving
     !> holds fields and parameters related to calving
     real(dp),dimension(:,:),  pointer :: calving_thck => null()   !> thickness loss in grid cell due to calving
                                                                   !> scaled by thk0 like mass balance, thickness, etc.
     real(dp),dimension(:,:),  pointer :: calving_rate => null()   !> rate of ice loss due to calving (m/yr ice)
     real(dp),dimension(:,:),  pointer :: calving_rate_tavg => null()  !> rate of ice loss due to calving (m/yr ice, time average)
     integer, dimension(:,:),  pointer :: calving_mask => null()   !> calve floating ice wherever the mask = 1 (whichcalving = CALVING_GRID_MASK)
     real(dp),dimension(:,:),  pointer :: lateral_rate => null()   !> lateral calving rate (m/yr, not scaled)
                                                                   !> (whichcalving = EIGENCALVING, CALVING_DAMAGE) 
     real(dp),dimension(:,:),  pointer :: tau_eigen1 => null()     !> first eigenvalue of 2D horizontal stress tensor (Pa)
     real(dp),dimension(:,:),  pointer :: tau_eigen2 => null()     !> second eigenvalue of 2D horizontal stress tensor (Pa)
     real(dp),dimension(:,:),  pointer :: tau_eff => null()        !> effective stress (Pa) for calving; derived from tau_eigen1, tau_eigen2
     real(dp),dimension(:,:,:),pointer :: damage => null()         !> 3D damage tracer, 0 > damage < 1 (whichcalving = CALVING_DAMAGE)
  
     real(dp) :: marine_limit =  -200.d0         !> value of topg/relx at which floating ice calves (m)
                                                 !> (whichcalving = CALVING_RELX_THRESHOLD, CALVING_TOPG_THRESHOLD)
     real(dp) :: calving_fraction = 0.2d0        !> fractional thickness of floating ice that calves
                                                 !> (whichcalving = CALVING_FLOAT_FRACTION)
                                                 !> WHL - previously defined as the fraction of floating ice that does not calve
     real(dp) :: timescale = 0.0d0               !> timescale (yr) for calving (Glissade only); calving_thck = thck*max(dt/calving_timescale,1)
                                                 !> if calving_timescale = 0, then the full column calves at once
     real(dp) :: minthck = 100.d0                !> minimum thickness (m) of floating ice at marine edge before it calves
                                                 !> (whichcalving = CALVING_THCK_THRESHOLD or EIGENCALVING)
     real(dp) :: eigencalving_constant = 0.01d0  !> eigencalving constant, lateral calving rate (m/yr) per unit stress (Pa)
                                                 !> (whichcalving = EIGENCALVING)
     real(dp) :: eigen2_weight = 1.0d0           !> weight given to tau_eigen2 relative to tau_eigen1 in tau_eff (unitless)
     real(dp) :: damage_constant = 1.0d-7        !> damage constant; rate of change of damage (1/yr) per unit stress (Pa)
                                                 !> (whichcalving = CALVING_DAMAGE) 
     real(dp) :: damage_threshold = 0.75d0       !> threshold at which ice column is deemed sufficiently damaged to calve
                                                 !> assuming that 0 = no damage, 1 = total damage (whichcalving = CALVING_DAMAGE)
     real(dp) :: lateral_rate_max = 3000.d0      !> max lateral calving rate (m/yr) for damaged ice (whichcalving = CALVING_DAMAGE)
     integer :: ncull_calving_front = 0          !> number of times to cull calving_front cells at initialization, if cull_calving_front = T
                                                 !> Set to a larger value to remove wider peninsulas
     real(dp) :: taumax_cliff = 1.0d6            !> yield stress (Pa) for marine-based ice cliffs
     real(dp) :: cliff_timescale = 10.0d0        !> time scale (yr) for limiting marine cliffs (yr)
     real(dp) :: calving_front_x = 0.0d0         !> for CALVING_GRID_MASK option, calve ice wherever abs(x) > calving_front_x (m)
     real(dp) :: calving_front_y = 0.0d0         !> for CALVING_GRID_MASK option, calve ice wherever abs(y) > calving_front_y (m)
                                                 !> NOTE: This option is applied only if calving_front_x or calving_front_y > 0

  end type glide_calving

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type eismint_climate_type

     ! holds parameters for the eismint climate

     ! For EISMINT2:
     ! airt(1) = Tmin = summit surface temperature (K)
     ! airt(2) = S_T = horizontal temperature gradient (K/m)
     ! nmsb(1) = M_max = max accumulation (m/yr)
     ! nmsb(2) = S_b = horizontal smb gradient (m/yr/m)
     ! nmsb(3) = R_el = radial distance from summit where mass balance = 0 (m)   
     !

     integer :: eismint_type = 0
     !> select EISMINT experiment
     !> \begin{description}
     !> \item[{\bf 1}] EISMINT-1 fixed margin
     !> \item[{\bf 2}] EISMINT-1 moving margin
     !> \item[{\bf 3}] EISMINT-2
     !> \item[{\bf 4}] MISMIP-1 (not EISMINT but has similar climate parameters)
     !> \item[{\bf 5}] Exact verification (not EISMINT but has similar climate parameters)
     !> \end{description}

     ! NOTE: The initial nmsb values in the declarations below are appropriate
     !       for EISMINT-2, but the initial airt values are not.
     ! TODO: Change default airt values in eismint_type to be consistent with EISMINT-2?

     !> air temperature parameterisation K, K km$^{-3}$
     real(dp), dimension(2) :: airt = (/ -3.15d0, 1.d-2 /)  

     !> mass balance parameterisation:
!!     real(dp), dimension(3) :: nmsb = (/ 0.5d0, 1.05d-5, 450.0d3 /)
     real(dp), dimension(3) :: nmsb = (/ 0.5d0, 1.0d-5, 450.0d3 /)

     !> EISMINT time-dep climate forcing period, switched off when set to 0
     real(dp) :: period = 0.d0

     !> EISMINT amplitude of mass balance time-dep climate forcing
     real(dp) :: mb_amplitude = 0.2d0

  end type eismint_climate_type

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_temper

    !> Holds fields relating to temperature.

    !Note: In the Glide dycore, temp, flwa and dissip live on the unstaggered vertical grid
    !       at layer interfaces and have vertical dimension (1:upn).
    !      In the Glissade dycore, with remapping advection of temperature, 
    !       temp, flwa and dissip live on the staggered vertical grid at layer midpoints.  
    !       The vertical dimensions are (0:upn) for temp and (1:upn-1) for flwa and dissip.
    !
    !      bheatflx, ucondflx, and lcondflx are defined as positive down,
    !       so they will often be < 0.  
    !      However, bfricflx and dissipcol are defined to be >= 0.
    !
    !      If bheatflx is read from a data file, be careful about the sign!
    !      In input data, the geothermal heat flux may be defined as positive upward,
    !       whereas bheatflx is defined as positive downward.

    real(dp),dimension(:,:,:),pointer :: temp => null()      !> 3D temperature field.
    real(dp),dimension(:,:,:),pointer :: tempunstag => null()!> 3D temperature field unstaggered in the vertical.
    real(dp),dimension(:,:),  pointer :: bheatflx => null()  !> basal heat flux (W/m^2) (geothermal, positive down)
    real(dp),dimension(:,:,:),pointer :: flwa => null()      !> Glen's flow factor $A$.
    real(dp),dimension(:,:,:),pointer :: dissip => null()    !> interior heat dissipation rate, divided by rhoi*Ci (deg/s)
    real(dp),dimension(:,:),  pointer :: btemp => null()     !> Basal temperature on ice grid; diagnosed from temp(upn)
    real(dp),dimension(:,:),  pointer :: stagbtemp => null() !> Basal temperature on velo grid
    real(dp),dimension(:,:),  pointer :: bpmp => null()      !> Basal pressure melting point temperature
    real(dp),dimension(:,:),  pointer :: stagbpmp => null()  !> Basal pressure melting point temperature on velo grid
    real(dp),dimension(:,:),  pointer :: bfricflx => null()  !> basal heat flux (W/m^2) from friction (>= 0)
    real(dp),dimension(:,:,:),pointer :: waterfrac => null() !> fractional water content in layer (0 <= waterfrac <= 1)
    real(dp),dimension(:,:,:),pointer :: enthalpy => null()  !> specific enthalpy in layer (J m-3)
                                                             !> = rhoi * Ci * T for cold ice
    !TODO - Remove ucondflx, lcondflx, dissipcol and make these local to glissade_therm?
    !       Probably cannot remove ucondflx because it may be needed for coupling.
    real(dp),dimension(:,:),  pointer :: ucondflx => null()  !> conductive heat flux (W/m^2) at upper sfc (positive down)
    real(dp),dimension(:,:),  pointer :: lcondflx => null()  !> conductive heat flux (W/m^2) at lower sfc (positive down)
    real(dp),dimension(:,:),  pointer :: dissipcol => null() !> total heat dissipation rate (W/m^2) in column (>= 0)

     ! fields related to basal water
     !TODO - Move these fields to the basal_physics type?
     real(dp),dimension(:,:),  pointer :: bwat => null()      !> Basal water depth
     real(dp),dimension(:,:),  pointer :: bwatflx => null()   !> Basal water flux 
     real(dp),dimension(:,:),  pointer :: stagbwat => null()  !> Basal water depth on velo grid

    real(dp) :: pmp_offset = 5.0d0        ! offset of initial Tbed from pressure melting point temperature (deg C)
    real(dp) :: pmp_threshold = 1.0d-3    ! bed is assumed thawed where Tbed >= pmptemp - pmp_threshold (deg C)

    !TODO - Remove some of the following from the derived type
    integer  :: niter   = 0   
    real(dp) :: perturb = 0.d0
    real(dp) :: grid    = 0.d0 
    integer  :: tpt     = 0      !> Pointer to time series data
    logical  :: first1  = .true. !>
    logical  :: newtemps = .false. !> new temperatures

  end type glide_temper

  type glide_inversion

     ! Notes on inversion fields:
     ! With which_ho_inversion = HO_INVERSION_COMPUTE, bmlt_float_inversion is computed and applied during each step.
     ! For which_ho_inversion = HO_INVERSION_PRESCRIBE, bmlt_float_prescribed (as computed in a
     !  previous inversion run) is read from the input file. 

     real(dp), dimension(:,:), pointer :: &
          bmlt_float_inversion => null(),      & !> basal melt rate computed by inversion;
                                                 !> used to relax thickness of floating ice toward observed target
          bmlt_float_prescribed => null()        !> basal melt rate prescribed from a previous inversion

     real(dp), dimension(:,:), pointer :: &
          powerlaw_c_inversion => null(),      & !> spatially varying powerlaw_c field, Pa (m/yr)^(-1/3)
          powerlaw_c_prescribed => null(),     & !> powerlaw_c field, prescribed from a previous inversion
          usrf_inversion => null(),            & !> upper surface elevation, used for Cp inversion (m)
          dthck_dt_inversion => null()           !> rate of thickness change, used for Cp inversion (m/s) 

     ! parameters for inversion of basal friction coefficients
     ! Note: These values work well for MISMIP+, but may not be optimal for whole ice sheets.
     ! Note: inversion_babc_timescale and inversion_babc_dthck_dt_scale are later rescaled to SI units (s and m/s).

     real(dp) ::  &
          powerlaw_c_max = 1.0d5,             &  !> max value of powerlaw_c, Pa (m/yr)^(-1/3)
          powerlaw_c_min = 1.0d2,             &  !> min value of powerlaw_c, Pa (m/yr)^(-1/3)
          powerlaw_c_land = 2.0d-4,           &  !> default value of powerlaw_c on land (topg >= eus)
          powerlaw_c_marine = 1.0d-3             !> default value of powerlaw_c below sea level

     real(dp) ::  &
          babc_timescale = 500.d0,        & !> inversion timescale (yr); must be > 0
          babc_thck_scale = 100.d0,       & !> thickness inversion scale (m); must be > 0
          babc_dthck_dt_scale = 0.10d0,   & !> dthck_dt inversion scale (m/yr); must be > 0
          babc_space_smoothing = 1.0d-2,  & !> factor for spatial smoothing of powerlaw_c; larger => more smoothing
          babc_time_smoothing = 0.0d0       !> factor for exponential moving average of usrf/dthck_dt_inversion
                                            !> range [0,1]; larger => slower discounting of old values, more smoothing

     real(dp) ::  &
          bmlt_thck_buffer = 1.0d0          !> cells with a grounded target are restored to thck_flotation + bmlt_thck_buffer (m)

  end type glide_inversion

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_basal_melt

     !> Holds fields and parameters relating to basal melting

     !Note: In the Glide dycore, the only active field in this type is bmlt.
     !      The other fields are used in Glissade only.

     !WHL - debug
     real(dp), dimension(:,:), pointer :: &
          bmlt_applied_old => null(),  &
          bmlt_applied_diff => null()

     ! bmlt fields for grounded and floating ice

     real(dp), dimension(:,:), pointer :: &
          bmlt => null(),                         & !> basal melt rate (> 0 for melt, < 0 for freeze-on)
                                                    !> bmlt = bmlt_ground + bmlt_float
          bmlt_applied => null(),                 & !> basal melt rate applied to ice (m/yr)
                                                    !>    = 0 for ice-free cells with bmlt > 0
          bmlt_applied_tavg => null(),            & !> basal melt rate applied to ice (m/yr, time average)
          bmlt_ground => null(),                  & !> basal melt rate for grounded ice
          bmlt_float => null(),                   & !> basal melt rate for floating ice
          bmlt_float_external => null(),          & !> external basal melt rate field
          bmlt_float_anomaly => null()              !> basal melt rate anomaly field

     real(dp) :: bmlt_float_factor = 1.0d0          !> adjustment factor for external bmlt_float field

     ! MISMIP+ parameters for Ice1 experiments (BMLT_FLOAT_MISMIP)
     ! Note: Parameters with units yr^{-1} are scaled to s^{-1} in subroutine glide_scale_params
     real(dp) :: bmlt_float_omega = 0.2d0           !> time scale for basal melting (yr-1)
                                                    !> default value = 0.2 yr^{-1} for MISMIP+ Ice1r
     real(dp) :: bmlt_float_h0 = 75.d0              !> scale for sub-shelf cavity thickness (m)
                                                    !> default value = 75 m for MISMIP+ Ice1r
     real(dp) :: bmlt_float_z0 = -100.d0            !> scale for ice draft, relative to sea level (m)
                                                    !> default value = -100 m for MISMIP+ Ice1r

     ! MISMIP+ parameters for Ice2 experiments (BMLT_FLOAT_CONSTANT)
     real(dp) :: bmlt_float_const = 0.d0            !> constant melt rate (m/yr)
                                                    !> set to 100 m/yr for MISMIP+ Ice2r
     real(dp) :: bmlt_float_xlim = 0.d0             !> melting is allowed only for abs(x1) > bmlt_float_xlim
                                                    !> set to 480 km for MISMIP+ Ice2r
     ! parameters for BMLT_FLOAT_DEPTH
     ! The sub-shelf melt rate is piecewise linear, generally with greater melting at depth.
     ! The maximum melting and freezing rates are set independently.
     ! The melting/freezing rates fall linearly from the max values to zero over ranges defined by
     !  zmeltmax, zmelt0 and zfrzmax.
     real(dp) :: bmlt_float_depth_meltmax = 10.d0     !> max melt rate at depth (m/yr)
     real(dp) :: bmlt_float_depth_frzmax = 0.d0       !> max freezing rate near surface (m/yr)
     real(dp) :: bmlt_float_depth_zmeltmax = -500.d0  !> depth (m) below which bmlt_float = meltmax
     real(dp) :: bmlt_float_depth_zmelt0 = -200.d0    !> depth (m) where bmlt_float = 0
     real(dp) :: bmlt_float_depth_zfrzmax = -100.d0   !> depth (m) above which bmlt_float = -frzmax

     ! parameters and fields for warm ocean
     ! Where warm_ocean_mask = 1, the sub-shelf melt rate is meltmin above zmeltmin, then increases linearly
     !  to meltmax at zmeltmax. 
     integer, dimension(:,:), pointer :: warm_ocean_mask => null()  !> mask for applying warm ocean melt rate
     real(dp) :: bmlt_float_depth_meltmin = 0.d0      !> min melt rate in warm ocean (m/yr)
     real(dp) :: bmlt_float_depth_zmeltmin = 0.d0     !> depth (m) above which bmlt_float = meltmin

     ! initMIP-Antarctica parameters
     real(dp) :: bmlt_anomaly_timescale = 0.0d0     !> number of years over which the bmlt_float anomaly is phased in linearly
                                                    !> If set to zero, then the anomaly is applied immediately.

  end type glide_basal_melt

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_plume

     !> Holds fields and parameters relating to a sub-shelf plume model
     !> Note: Entrainment/detrainment rates are computed with units of m/s but output with m/yr

     real(dp),dimension(:,:), pointer :: T_basal => null()       !> basal ice temperature; at freezing point (deg C)
     real(dp),dimension(:,:), pointer :: S_basal => null()       !> basal salinity; at freezing point (psu)
     real(dp),dimension(:,:), pointer :: u_plume => null()       !> x component of plume velocity at cell centers (m/s)
     real(dp),dimension(:,:), pointer :: v_plume => null()       !> y component of plume velocity at cell centers (m/s)
                                                                 !> Note: Plume velocities are prognosed on edges, then interpolated
                                                                 !>       to cell centers for diagnostics 
     real(dp),dimension(:,:), pointer :: u_plume_Cgrid => null() !> x component of plume velocity on C grid edges (m/s)
     real(dp),dimension(:,:), pointer :: v_plume_Cgrid => null() !> y component of plume velocity on C grid edges (m/s)
     real(dp),dimension(:,:), pointer :: D_plume => null()       !> plume thickness (m)
     real(dp),dimension(:,:), pointer :: ustar_plume => null()   !> plume friction velocity (m/s) on ice grid
     real(dp),dimension(:,:), pointer :: drho_plume => null()    !> density difference between plume and ambient ocean (kg/m3)
     real(dp),dimension(:,:), pointer :: T_plume => null()       !> plume temperature (deg C)
     real(dp),dimension(:,:), pointer :: S_plume => null()       !> plume salinity (psu)
     real(dp),dimension(:,:), pointer :: entrainment => null()   !> entrainment rate from ambient ocean to plume (positive up)
     real(dp),dimension(:,:), pointer :: detrainment => null()   !> detrainment rate from plume to ambient ocean (positive down)
     real(dp),dimension(:,:), pointer :: divDu_plume => null()   !> divergence of D_plume*u_plume
     real(dp),dimension(:,:), pointer :: T_ambient => null()     !> ambient ocean temperature below ice and plume (deg C)
     real(dp),dimension(:,:), pointer :: S_ambient => null()     !> ambient ocean salinity below ice and plume (psu)

     ! MISOMIP parameters
     ! Note: T0, Tbot, S0, Sbot, gammaT and gammaS can be set in the config file.
     !     - T0, Tbot, S0 and Sbot are prescribed for cold and warm profiles.
     !     - gammaT is to be tuned to give a bmlt_float of the desired mean value.
     !     - gammaS should equal gammaS/35.
     real(dp) :: T0 = -1.9d0             !> sea surface temperature (deg C)
     real(dp) :: Tbot =  1.0d0           !> temperature at the sea floor, warm profile (deg C)
!     real(dp) :: Tbot = -1.9d0          !> temperature at the sea floor, cold profile (deg C)
     real(dp) :: S0 = 33.8d0             !> sea surface salinity (psu)
     real(dp) :: Sbot = 34.7d0           !> salinity at the sea floor, warm profile (psu)
!     real(dp) :: Sbot = 34.55d0         !> salinity at the sea floor, cold profile (psu)
     real(dp) :: gammaT = 2.2d-2         !> nondimensional heat transfer coefficient
                                         !> value of 2.2e-2 suggested by Asay-Davis et al. as an initial guess
     real(dp) :: gammaS = 2.2d-2/35.d0   !> nondimensional salt transfer coefficient
                                         !> for MISOMIP, should be set to gammaT/35 
     real(dp) :: zbed_deep = -720.d0     !> min sea floor elevation (m)

  end type glide_plume

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_basal_physics
     !> Holds variables related to basal physics associated with ice dynamics
     !> See glissade_basal_traction.F90 for usage details

     !Note: By default, beta_grounded_min is set to a small nonzero value.
     !      Larger values (~10 to 100 Pa yr/m) might be needed for stability in realistic simulations.
     real(dp) :: beta_grounded_min = 1.0d0  !> minimum value of beta for grounded ice, Pa yr/m (glissade only; scaled during init)
     
     real(dp) :: ho_beta_const = 1000.d0    !> spatially uniform beta for HO dycores, Pa yr/m (scaled during init)
     real(dp) :: ho_beta_small = 1000.d0    !> small beta for sliding over a thawed bed, Pa yr/m (scaled during init)
     real(dp) :: ho_beta_large = 1.0d10     !> large beta to enforce (virtually) no slip, Pa yr/m (scaled during init)

     integer,  dimension(:,:), pointer :: bpmp_mask => null()   !> basal pressure melting point mask; = 1 where Tbed = bpmp, elsewhere = 0
                                                                !> Note: Defined on velocity grid, whereas temp and bpmp are on ice grid

     ! Note: It may make sense to move effecpress to a hydrology model when one is available.
     real(dp), dimension(:,:), pointer :: effecpress => null()          !> effective pressure (Pa)
     real(dp), dimension(:,:), pointer :: effecpress_stag => null()     !> effective pressure on staggered grid (Pa)
     real(dp), dimension(:,:), pointer :: C_space_factor => null()      !> spatial factor for basal shear stress (no dimension)
     real(dp), dimension(:,:), pointer :: C_space_factor_stag => null() !> spatial factor for basal shear stress on staggered grid (no dimension)
     real(dp), dimension(:,:), pointer :: tau_c => null()               !> yield stress for plastic sliding (Pa)

     ! parameters for reducing the effective pressure where the bed is warm, saturated or connected to the ocean
     real(dp) :: effecpress_delta = 0.02d0          !> multiplier for effective pressure N where the bed is saturated and/or thawed (unitless)
     real(dp) :: effecpress_bpmp_threshold = 0.1d0  !> temperature range over which N ramps from a small value to full overburden (deg C)
     real(dp) :: effecpress_bmlt_threshold = 1.0d-3 !> basal melting range over which N ramps from a small value to full overburden (m/yr)
     real(dp) :: p_ocean_penetration = 0.0d0        !> p-exponent parameter for ocean penetration parameterization (unitless, 0 <= p <= 1)

     ! parameters for pseudo-plastic sliding law (based on PISM)
     ! (tau_bx,tau_by) = -tau_c * (u,v) / (u_0^q * |u|^(1-q))
     ! where the yield stress tau_c = tan(phi) * N
     ! N = effective pressure

     real(dp) :: pseudo_plastic_q = 0.5d0        !> exponent for pseudo-plastic law (unitless), 0 <= q <= 1
                                                 !> q = 1 => linear sliding law; q = 0 => plastic; intermediate values => power law
     real(dp) :: pseudo_plastic_u0 = 100.d0      !> threshold velocity for pseudo-plastic law (m/yr)

     ! The following 4 parameters give a linear increase in phi between elevations bedmin and bedmax
     real(dp) :: pseudo_plastic_phimin =    5.d0 !> min(phi) in pseudo-plastic law, for topg <= bedmin (degrees, 0 < phi < 90)
     real(dp) :: pseudo_plastic_phimax =   40.d0 !> max(phi) in pseudo-plastic law, for topg >= bedmax (degrees, 0 < phi < 90)
     real(dp) :: pseudo_plastic_bedmin = -700.d0 !> bed elevation (m) below which phi = phimin
     real(dp) :: pseudo_plastic_bedmax =  700.d0 !> bed elevation (m) above which phi = phimax

     ! parameters for friction powerlaw
     real(dp) :: friction_powerlaw_k = 8.4d-9    !> coefficient (m y^-1 Pa^-2) for the friction power law based on effective pressure
                                                 !> The default value is from Bindschadler (1983) based on fits to observations, converted to CISM units.

     ! parameters for Coulomb friction sliding law (default values from Pimentel et al. 2010)
     real(dp) :: coulomb_c = 0.42d0              !> basal stress constant (no dimension)
                                                 !> Pimentel et al. have coulomb_c = 0.84*m_max, where m_max = coulomb_bump_max_slope
     real(dp) :: coulomb_bump_wavelength = 2.0d0 !> bedrock wavelength at subgrid scale precision (m)
     real(dp) :: coulomb_bump_max_slope = 0.5d0  !> maximum bed bump slope at subgrid scale precision (no dimension)
     real(dp) :: flwa_basal = 1.0d-16            !> Glen's A at the bed for Schoof (2005) Coulomb friction law (Pa^{-n} yr^{-1})
                                                 !> = 3.1688d-24 Pa{-n} s{-1}, the value used by Leguy et al. (2014)

     ! parameters for power law, taub_b = C * u_b^(1/m); used for HO_BABC_COULOMB_POWERLAW_TSAI/SCHOOF
     ! The default values are from Asay-Davis et al. (2016).
     ! The value of powerlaw_c suggested by Tsai et al. (2015) is 7.624d6 Pa m^(-1/3) s^(1/3).
     ! This value can be converted to CISM units by dividing by scyr^(1/3), to obtain 2.413d4 Pa m^(-1/3) yr^(1/3).
     ! Note: The Tsai et al. Coulomb friction law uses coulomb_c above, with
     !       effective pressure N as in Leguy et al. (2014) with p_ocean_penetration = 1.

     real(dp) :: powerlaw_c = 1.0d4              !> friction coefficient in power law, units of Pa m^(-1/3) yr^(1/3)
     real(dp) :: powerlaw_m = 3.d0               !> exponent in power law (unitless)
      
     ! parameter to limit the min value of beta for various power laws
     real(dp) :: beta_powerlaw_umax = 0.0d0      !> upper limit of ice speed (m/yr) when evaluating powerlaw beta
                                                 !> Where u > umax, let u = umax when evaluating beta(u)

     ! parameter for constant basal water
     ! Note: This parameter applies to HO_BWAT_CONSTANT only.
     !       For Glide's BWATER_CONST, the constant value is hardwired in subroutine calcbwat.
     real(dp) :: const_bwat = 10.d0              !> constant basal water depth (m)

     ! parameters for local till model
     ! The default values are from Aschwanden et al. (2016) and Bueler and van Pelt (2015).
     real(dp) :: bwat_till_max = 2.0d0           !> maximum water depth in till (m)
     real(dp) :: c_drainage = 1.0d-3             !> uniform drainage rate (m/yr)
     real(dp) :: N_0 = 1000.d0                   !> reference effective pressure (Pa)
     real(dp) :: e_0 = 0.69d0                    !> reference void ratio (dimensionless)
     real(dp) :: C_c = 0.12d0                    !> till compressibility (dimensionless)
                                                 !> Note: The ratio (e_0/C_c) is the key parameter

     ! Note: A basal process model is not currently supported, but a specified mintauf can be passed to subroutine calcbeta
     !       to simulate a plastic bed..
     real(dp),dimension(:,:)  ,pointer :: mintauf => null() ! Bed strength (yield stress) calculated with basal process model

  end type glide_basal_physics

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_lithot_type
     !> holds variables for temperature calculations in the lithosphere

     real(dp),dimension(:,:,:),pointer :: temp => null()    !> Three-dimensional temperature field.
     logical, dimension(:,:), pointer :: mask => null()     !> whether the point has been ice covered at some time

     integer :: num_dim = 1                                 !> either 1 or 3 for 1D/3D calculations

     ! The sparse matrix and linearised arrays
     type(sparse_matrix_type) :: fd_coeff, fd_coeff_slap
     integer :: all_bar_top
     real(dp), dimension(:), pointer :: rhs
     real(dp), dimension(:), pointer :: answer
     real(dp), dimension(:), pointer :: supd,diag,subd

     ! work arrays for solver
     real(dp), dimension(:), pointer :: rwork
     integer, dimension(:), pointer :: iwork
     integer mxnelt

     real(dp), dimension(:), pointer :: deltaz => null()    !> array holding grid spacing in z
     real(dp), dimension(:,:), pointer :: zfactors => null()!> array holding factors for finite differences of vertical diffu
     real(dp) :: xfactor,yfactor !> factors for finite differences of horizontal diffu


     real(dp) :: surft = 2.d0          !> surface temperature, used for calculating initial temperature distribution
     real(dp) :: mart  = 2.d0          !> sea floor temperature 
     integer  :: nlayer = 20           !> number of layers in lithosphere
     real(dp) :: rock_base = -5000.d0  !> depth below sea-level at which geothermal heat gradient is applied
     
     integer :: numt = 0        !> number time steps for spinning up GTHF calculations

     real(dp) :: rho_r = 3300.0d0 !> The density of lithosphere (kg m$^{-3}$)
     real(dp) :: shc_r = 1000.0d0 !> specific heat capcity of lithosphere (J kg$^{-1}$ K$^{-1}$)
     real(dp) :: con_r = 3.3d0    !> thermal conductivity of lithosphere (W m$^{-1}$ K$^{-1}$)

     real(dp) :: diffu = 0. !> diffusion coefficient

  end type glide_lithot_type

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type isos_elastic
     !> Holds data used by isostatic adjustment calculations

     real(dp) :: d = 0.24d25                !> flexural rigidity (N m)
     real(dp) :: lr                         !> radius of relative stiffness
     real(dp) :: a                          !> radius of disk
     real(dp) :: c1,c2,cd3,cd4              !> coefficients
     real(dp), dimension(:,:), pointer :: w !> matrix operator for lithosphere deformation
     integer :: wsize                       !> size of operator (0:rbel_wsize, 0:rbel_wsize), operator is axis symmetric
  end type isos_elastic

  type isostasy_type
     !> contains isostasy configuration

     integer :: lithosphere = 1
     !> method for calculating equilibrium bedrock depression
     !> \begin{description}
     !> \item[0] local lithosphere, equilibrium bedrock depression is found using Archimedes' principle
     !> \item[1] elastic lithosphere, flexural rigidity is taken into account
     !> \end{description}

     integer :: asthenosphere = 1
     !> method for approximating the mantle
     !> \begin{description}
     !> \item[0] fluid mantle, isostatic adjustment happens instantaneously
     !> \item[1] relaxing mantle, exponential adjustment toward (relx - load)
     !> \end{description}

    integer :: whichrelaxed = 0

    !> relaxed topography:
    !> \begin{description}
    !> \item[0] get current topo (topg) and relaxed topo (relx) from separate input fields
    !> \item[1] first time slice of input topo is relaxed
    !> \item[2] first time slice of input topo is in isostatic equilibrium
    !> \end{description}

     real(dp) :: relaxed_tau = 4000.d0        ! characteristic time constant of relaxing mantle (yr)
     real(dp) :: period = 100.d0              ! lithosphere update period (yr)
     integer :: nlith                         ! update lithosphere every nlith time steps; tinc * nlith = lithosphere_period
     logical :: new_load = .false.            ! set to true if there is a new surface load
     type(isos_elastic) :: rbel               ! structure holding elastic lithosphere setup

     real(dp),dimension(:,:),pointer :: relx => null()  ! elevation of relaxed topography, m/thk0
     real(dp),dimension(:,:),pointer :: load => null()  ! deflection due to applied load on lithosphere, m/thk0
                                                        ! defined as positive for downward deflection
     real(dp),dimension(:,:),pointer :: load_factors => null() ! temporary used for load calculation

  end type isostasy_type

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_funits
    character(fname_length) :: sigfile=''                      !> sigma coordinates file
    character(fname_length) :: ncfile=''                       !> configuration file for netCDF I/O
    type(glimmer_nc_output),pointer :: out_first=>NULL()       !> first element of linked list defining netCDF outputs
    type(glimmer_nc_input), pointer :: in_first=>NULL()        !> first element of linked list defining netCDF inputs
    type(glimmer_nc_input), pointer :: frc_first=>NULL()       !> first element of linked list defining netCDF forcings
    ! Note: forcing files are of the same type as input files since they share a lot in common.
  end type glide_funits

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_numerics

    !> Parameters relating to the model numerics
    real(dp) :: tstart =     0.d0 !> starting time
                                  !> Note: For restarts, tstart is set to the time the restart begins
    real(dp) :: tend   =  1000.d0 !> end time
    real(dp) :: time   =     0.d0 !> main loop counter in years
    real(dp) :: tinc   =     1.d0 !> time step of main loop in years 
    real(dp) :: ntem   =     1.d0 !> multiplier of main time step; allows longer temperature time step
    real(dp) :: alpha  =    0.5d0 !> richard suggests 1.5 - was a parameter in original
    real(dp) :: alphas =    0.5d0 !> was a parameter in the original
    real(dp) :: thklim =   100.d0 !> min thickness for computing ice dynamics (m) 
    real(dp) :: thklim_temp = 1.d0    !> min thickness for computing vertical temperature (m) (higher-order only)
    real(dp) :: thck_gradient_ramp = 0.d0 !> thickness scale over which gradients increase from zero to full value (HO only)
    real(dp) :: dew    =    20.d3     !> grid cell size in east-west direction
    real(dp) :: dns    =    20.d3     !> grid cell size in north-south direction
    real(dp) :: dt     =     0.d0     !> ice dynamics timestep
    real(dp) :: dttem  =     0.d0     !> temperature timestep
    real(dp) :: dt_transport = 0.d0   !> timestep for subcycling transport within the dynamics timestep dt
    real(dp) :: nshlf  =     0.d0          !TODO - not currently used; remove?
    integer  :: subcyc =     1
    real(dp) :: periodic_offset_ew = 0.d0 !> optional periodic_offsets for ismip-hom and similar tests;
    real(dp) :: periodic_offset_ns = 0.d0 !> may be needed to ensure continuous ice geometry at edges of the global domain
    integer  :: tstep_count = 0   !> number of time steps since the start of the simulation
    
    ! Vertical coordinate ---------------------------------------------------
                                                               
    real(dp),dimension(:),pointer :: sigma => null() !> Sigma values for vertical spacing of 
                                                     !> model levels
    real(dp),dimension(:),pointer :: stagsigma => null() !> Staggered values of sigma (layer midpts)
    real(dp),dimension(:),pointer :: stagwbndsigma => null() !> Staggered values of sigma (layer midpts) with boundaries

    integer :: profile_period = 100            ! profile frequency

    real(dp) :: dt_diag = 0.d0     ! diagnostic interval (write diagnostic output every dt_diag years)
                                   ! dt_diag = 0 => never write diagnostic output
    integer  :: ndiag = 0          ! diagnostic interval (write diagnostic output every ndiag timesteps)
                                   ! ndiag = 0 => never write diagnostic output
    integer  :: idiag = 1          ! global grid indices for diagnostic point
    integer  :: jdiag = 1          ! 
    integer  :: idiag_local = 1    ! local grid indices for diagnostic point
    integer  :: jdiag_local = 1
    integer  :: rdiag_local = 0    ! task number for diagnostic point

    ! CFL diagnostics and adaptive subcycling
    real(dp) :: adv_cfl_dt = 0.0d0  ! maximum allowable dt (yrs) based on advective CFL (calculated by model for each time step)
    real(dp) :: diff_cfl_dt = 0.0d0 ! maximum allowable dt (yrs) based on diffusive CFL (calculated by model for each time step)
    real(dp) :: adaptive_cfl_threshold = 0.0d0  ! adaptively subcycle the advection when advective CFL exceeds this value
                                                ! (zero value means no adaptive subcycling)

  end type glide_numerics

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_velowk
    real(dp),dimension(:),  pointer :: depth    => null()
    real(dp),dimension(:),  pointer :: dupsw    => null()
    real(dp),dimension(:),  pointer :: depthw   => null()
    real(dp),dimension(:),  pointer :: suvel    => null()
    real(dp),dimension(:),  pointer :: svvel    => null()
    real(dp),dimension(:,:),pointer :: fslip    => null()
    real(dp),dimension(:,:),pointer :: dintflwa => null()
    real(dp),dimension(:),  pointer :: dups     => null()
    real(dp),dimension(4) :: fact
    real(dp),dimension(4) :: c = 0.d0
    real(dp) :: watwd  = 3.0d0
    real(dp) :: watct  = 10.0d0
    real(dp) :: trc0   = 0.d0
    real(dp) :: trcmin = 0.0d0
    real(dp) :: marine = 1.0d0
    real(dp) :: trcmax = 10.0d0
    real(dp) :: btrac_const = 0.0d0  !TODO - Remove from glide_velowk type; already in glide_paramets type.
    real(dp) :: btrac_slope = 0.0d0
    real(dp) :: btrac_max = 0.d0
  end type glide_velowk

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_thckwk
     real(dp),dimension(:,:),  pointer :: oldthck   => null()
     real(dp),dimension(:,:),  pointer :: oldthck2  => null()
     real(dp),dimension(:,:,:),pointer :: olds      => null()
     integer  :: nwhich  = 2
     real(dp) :: oldtime = 0.d0
     
     ! next four are for ADI evolution only
     real(dp), dimension(:), pointer :: alpha => null()
     real(dp), dimension(:), pointer :: beta  => null()
     real(dp), dimension(:), pointer :: gamma => null()
     real(dp), dimension(:), pointer :: delta => null()

  end type glide_thckwk

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_tempwk
    real(dp),dimension(:,:,:),pointer :: inittemp => null()
    real(dp),dimension(:,:,:),pointer :: compheat => null()
    real(dp),dimension(:,:,:),pointer :: initadvt => null()
    real(dp),dimension(:),    pointer :: dupa     => null()
    real(dp),dimension(:),    pointer :: dupb     => null()
    real(dp),dimension(:),    pointer :: dupc     => null()
    real(dp),dimension(:),    pointer :: c1       => null()
    real(dp),dimension(:,:),  pointer :: dups     => null()
    real(dp),dimension(:,:),  pointer :: wphi     => null()
    real(dp),dimension(:,:),  pointer :: smth     => null()
    real(dp),dimension(:,:,:),pointer :: hadv_u   => null()
    real(dp),dimension(:,:,:),pointer :: hadv_v   => null()

    !*sfp** added space to the next 2 (cons, f) for use w/ HO and SSA dissip. calc. 
    real(dp),dimension(5)             :: cons     = 0.d0
    real(dp),dimension(5)             :: f        = 0.d0
    real(dp),dimension(8)             :: c        = 0.d0
    real(dp),dimension(2)             :: slide_f
    real(dp) :: noflow      = -1
    real(dp),dimension(2) :: advconst = 0.d0
    real(dp) :: zbed        = 0.d0
    real(dp) :: dupn        = 0.d0
    real(dp) :: wmax        = 0.d0
    real(dp) :: dt_wat      = 0.d0
    real(dp) :: watvel      = 0.d0
    integer  :: nwat        = 0
  end type glide_tempwk

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_paramets
    real(dp),dimension(5) :: bpar = (/ 0.2d0, 0.5d0, 0.0d0 ,1.0d-2, 1.0d0/)
    real(dp) :: btrac_const = 0.d0     ! m yr^{-1} Pa^{-1} (gets scaled during init)
    real(dp) :: btrac_slope = 0.0d0    ! Pa^{-1} (gets scaled during init)
    real(dp) :: btrac_max = 0.d0       ! m yr^{-1} Pa^{-1} (gets scaled during init)
    real(dp) :: geot   = -5.0d-2       ! W m^{-2}, positive down
    real(dp) :: flow_enhancement_factor = 1.0d0   ! flow enhancement parameter for the Arrhenius relationship;
                                                  ! typically > 1 for SIA models to speed up the ice
                                                  ! (Note the change relative to prev. versions of code - used to be 3.0)
    real(dp) :: flow_enhancement_factor_float = 1.0d0 ! flow enhancement parameter for floating ice
                                                      ! Default is 1.0, but for marine simulations a smaller value
                                                      !  may be needed to match observed shelf speeds
    real(dp) :: slip_ratio = 1.0d0     ! Slip ratio, used only in higher order code when the slip ratio beta computation is requested
    real(dp) :: hydtim = 1000.0d0      ! years, converted to s^{-1} and scaled
                                       ! 0 if no drainage
    real(dp) :: bwat_smooth = 0.01d0   ! basal water field smoothing strength

    !TODO - Move the following higher-order velocity parameters to a different type?
    real(dp) :: default_flwa = 1.0d-16 ! Glen's A to use in isothermal case, in units Pa^{-n} yr^{-1} 
                                       ! (would change to e.g. 4.6e-18 in EISMINT-ROSS case)
    real(dp) :: efvs_constant = 2336041.d0  ! value of efvs to use in constant efvs case, in units Pa yr
                                       ! = 0.5*A^(-1), where A = 2.140373 Pa^(-1) yr^(1) is the value used in ISMIP-HOM Test F
    real(dp) :: max_slope = 1.0d0      ! maximum surface slope allowed in Glissade dycore (unitless)
                                       ! Note: It may be necessary to reduce max_slope to ~0.1 to prevent huge velocities
                                       !       in regions of rough coastal topography
            
  end type glide_paramets

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  !TODO - Remove the basalproc type?  Not currently used.
!!  type glide_basalproc
    !Tuneables, set in the config file 
!    real(dp):: fric=0.45d0                   ! Till coeff of internal friction: ND
!    real(dp):: etillo=0.7d0                  ! Till void ratio at No
!    real(dp):: No=1000.d0                    ! Reference value of till effective stress
!    real(dp):: Comp=0.12d0                   ! Till coeff of compressibility: ND
!    real(dp):: Cv = 1.0d-8                   ! Till hydraulic diffusivity: m2/s
!    real(dp):: Kh = 1.0d-10                  !Till hydraulic conductivity: m/s
!    real(dp):: Zs = 3.0d0                    ! Solid till thickness: m
!    real(dp):: aconst=994000000d0            ! Constant in till strength eq. (Pa)
!    real(dp):: bconst=21.7d0                 ! Constant in till strength eq. (ND)
!    integer:: till_hot = 0
!    integer:: tnodes = 5

!    real(dp), dimension (:) , pointer :: till_dz => null()  !holds inital till layer spacing - 
    
    !Model variables that will be passed to other subroutines
!    real(dp),dimension(:,:)  ,pointer :: Hwater  => null() !Water available from till layer (m)
    !Model variables necessary for restart
!    real(dp),dimension(:,:,:)  ,pointer :: u => null()     !Till excess pore pressure (Pa)
!    real(dp),dimension(:,:,:)  ,pointer :: etill  => null()  !Till void ratio (ND)  
    
!  end type glide_basalproc

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  type glide_prof_type
     integer :: geomderv
     integer :: hvelos
     integer :: ice_mask1
     integer :: temperature
     integer :: ice_evo
     integer :: ice_mask2
     integer :: isos_water
     integer :: isos
  end type glide_prof_type

  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!TODO - Remove the glide_phaml type?  Commented out for now
!!  type glide_phaml
!!    real(dp),dimension(:,:),pointer :: uphaml => null()
!!    real(dp),dimension(:,:),pointer :: init_phaml => null()
!!    real(dp),dimension(:,:),pointer :: rs_phaml => null()
!!    !maybe put the x/y vectors here too just for simplicity
!!  end type glide_phaml
  
  !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  ! for JFNK, NOX Trilinos solver
  type, public :: glissade_solver

    integer ,dimension(:,:) ,allocatable :: ui 
    integer ,dimension(:,:) ,allocatable :: um 
    real(dp) ,dimension(:,:) ,allocatable :: d2thckcross
    real(dp) ,dimension(:,:) ,allocatable :: d2usrfcross
    integer ,dimension(2) :: pcgsize
    integer ,dimension(:) ,allocatable :: gxf
    real(dp)  :: L2norm
    type(sparse_matrix_type)  :: matrix
    type(sparse_matrix_type)  :: matrixA
    type(sparse_matrix_type)  :: matrixC
    real(dp),dimension(:),pointer :: rhsd    => null()
    real(dp),dimension(:),pointer :: answ    => null()
    integer :: ct     = 0

    !TODO - KJE - Remove ewn/nsn from glissade_solver type once new glide_global_type is working and we can use those ewn/nsn
    integer :: ewn
    integer :: nsn

  end type glissade_solver

       
  type glide_global_type    ! type containing all of the above for an ice sheet model instance
    integer              :: model_id !> Used in the global model list for error handling purposes
    type(glide_general)  :: general
    type(glide_options)  :: options
    type(glide_geometry) :: geometry
    type(glide_geomderv) :: geomderv
    type(glide_velocity) :: velocity
    type(glide_stress_t) :: stress   
    type(glide_climate)  :: climate
    type(eismint_climate_type) :: eismint_climate
    type(glide_calving)  :: calving
    type(glide_temper)   :: temper
    type(glide_basal_physics):: basal_physics
    type(glide_basal_melt)   :: basal_melt
    type(glide_inversion):: inversion
    type(glide_plume)    :: plume
    type(glide_lithot_type)  :: lithot
    type(glide_funits)   :: funits
    type(glide_numerics) :: numerics
    type(glide_velowk)   :: velowk
    type(glide_thckwk)   :: thckwk
    type(glide_tempwk)   :: tempwk
    type(glide_paramets) :: paramets
    type(glimmap_proj)   :: projection
    type(profile_type)   :: profile
    type(glide_prof_type):: glide_prof
    type(isostasy_type)  :: isostasy
    type(glissade_solver):: solver_data
!!    type(glide_basalproc):: basalproc
!!    type(glide_phaml)    :: phaml

  end type glide_global_type

contains

  subroutine glide_allocarr(model)    

    !> Allocates the model arrays, and initialises some of them to zero.
    !> These are the arrays allocated, and their dimensions:

    !TODO - Make sure the itemized lists in subroutine glide_allocarr are complete.

    !> In \texttt{model\%temper}:
    !> \begin{itemize}
    !> \item \texttt{temp(upn,0:ewn+1,0:nsn+1))}   !WHL - 2 choices
    !> \item \texttt{bheatflx(ewn,nsn))}
    !> \item \texttt{flwa(upn,ewn,nsn))}           !WHL - 2 choices
    !> \item \texttt{dissip(upn,ewn,nsn))}         !WHL - 2 choices
    !> \item \texttt{bwat(ewn,nsn))}
    !> \item \texttt{bfricflx(ewn,nsn))}
    !> \item \texttt{ucondflx(ewn,nsn))}
    !> \item \texttt{lcondflx(ewn,nsn))}
    !> \item \texttt{dissipcol(ewn,nsn))}
    !> \item \texttt{waterfrac(upn-1,ewn,nsn))}
    !> \item \texttt{enthalpy(0:upn,ewn,nsn))}
    !> \end{itemize}

    !> In \texttt{model\%basal_melt}:
    !> \begin{itemize}
    !> \item \texttt{bmlt(ewn,nsn))}
    !> \item \texttt{bmlt_ground(ewn,nsn))}
    !> \item \texttt{bmlt_applied(ewn,nsn))}
    !> \item \texttt{bmlt_float(ewn,nsn)}
    !> \item \texttt{bmlt_float_external(ewn,nsn)}
    !> \item \texttt{bmlt_float_anomaly(ewn,nsn)}
    !> \end{itemize}

    !> In \texttt{model\%inversion}:
    !> \item \texttt{bmlt_float_inversion(ewn,nsn)}
    !> \item \texttt{bmlt_float_prescribed(ewn,nsn)}
    !> \item \texttt{powerlaw_c_inversion(ewn,nsn)}
    !> \item \texttt{powerlaw_c_prescribed(ewn,nsn)}
    !> \item \texttt{usrf_inversion(ewn,nsn)}
    !> \item \texttt{dthck_dt_inversion(ewn,nsn)}

    !> In \texttt{model\%plume}:
    !> \begin{itemize}
    !> \item \texttt{T_basal(ewn,nsn)}
    !> \item \texttt{S_basal(ewn,nsn)}
    !> \item \texttt{u_plume(ewn,nsn)}
    !> \item \texttt{v_plume(ewn,nsn)}
    !> \item \texttt{u_plume_Cgrid(ewn,nsn)}
    !> \item \texttt{v_plume_Cgrid(ewn,nsn)}
    !> \item \texttt{D_plume(ewn,nsn)}
    !> \item \texttt{ustar_plume(ewn,nsn)}
    !> \item \texttt{drho_plume(ewn,nsn)}
    !> \item \texttt{T_plume(ewn,nsn)}
    !> \item \texttt{S_plume(ewn,nsn)}
    !> \item \texttt{entrainment(ewn,nsn)}
    !> \item \texttt{detrainment(ewn,nsn)}
    !> \item \texttt{divDu_plume(ewn,nsn)}
    !> \item \texttt{T_ambient(ewn,nsn)}
    !> \item \texttt{S_ambient(ewn,nsn)}
    !> \end{itemize}

    !> In \texttt{model\%velocity}:
    !> \begin{itemize}
    !> \item \texttt{uvel(upn,ewn-1,nsn-1))}
    !> \item \texttt{vvel(upn,ewn-1,nsn-1))}
    !> \item \texttt{wvel(upn,ewn,nsn))}
    !> \item \texttt{wgrd(upn,ewn,nsn))}
    !> \item \texttt{uflx(ewn-1,nsn-1))}
    !> \item \texttt{vflx(ewn-1,nsn-1))}
    !> \item \texttt{diffu(ewn,nsn))}
    !> \item \texttt{btrc(ewn,nsn))}
    !> \item \texttt{ubas(ewn,nsn))}
    !> \item \texttt{vbas(ewn,nsn))}
    !> \end{itemize}

    !> In \texttt{model\%climate}:
    !> \begin{itemize}
    !> \item \texttt{acab(ewn,nsn))}
    !> \item \texttt{artm(ewn,nsn))}
    !> \end{itemize}

    !> In \texttt{model\%geomderv}:
    !> \begin{itemize}
    !> \item \texttt{dthckdew(ewn,nsn))}
    !> \item \texttt{dusrfdew(ewn,nsn))}
    !> \item \texttt{dthckdns(ewn,nsn))}
    !> \item \texttt{dusrfdns(ewn,nsn))}
    !> \item \texttt{dthckdtm(ewn,nsn))}
    !> \item \texttt{dusrfdtm(ewn,nsn))}
    !> \item \texttt{stagthck(ewn-1,nsn-1))}
    !> \end{itemize}
  
    !> In \texttt{model\%geometry}:
    !> \begin{itemize}
    !> \item \texttt{thck(ewn,nsn))}
    !> \item \texttt{usrf(ewn,nsn))}
    !> \item \texttt{lsrf(ewn,nsn))}
    !> \item \texttt{topg(ewn,nsn))}
    !> \item \texttt{thck_obs(ewn,nsn))}
    !> \item \texttt{usrf_obs(ewn,nsn))}
    !> \item \texttt{topg_obs(ewn,nsn))}
    !> \item \texttt{mask(ewn,nsn))}
    !> \item \texttt{age(upn-1,ewn,nsn))}
    !> \item \texttt{tracers(ewn,nsn,ntracers,upn-1)}
    !> \item \texttt{f_flotation(ewn,nsn)}
    !> \item \texttt{f_ground(ewn-1,nsn-1)}
    !> \item \texttt{f_ground_cell(ewn-1,nsn-1)}
    !* (DFM) added floating_mask, ice_mask, lower_cell_loc, and lower_cell_temp
    !> \item \texttt{ice_mask(ewn,nsn))}
    !> \item \texttt{floating_mask(ewn,nsn))}
    !> \item \texttt{grounded_mask(ewn,nsn))}
    !> \item \texttt{ice_sheet_mask(ewn,nsn))}
    !> \item \texttt{ice_cap_mask(ewn,nsn))}
    !> \item \texttt{lower_cell_loc(ewn,nsn))}
    !> \item \texttt{lower_cell_temp(ewn,nsn))}
    !> \end{itemize}

    !> In \texttt{model\%thckwk}:
    !> \begin{itemize}
    !> \item \texttt{olds(ewn,nsn,thckwk\%nwhich))}
    !> \end{itemize}

    !> In \texttt{model\%numerics}:
    !> \begin{itemize}
    !> \item \texttt{sigma(upn))}
    !> \end{itemize}

    !> In \texttt{model\%numerics}:
    !> \begin{itemize}
    !> \item \texttt{stagsigma(upn-1))}
    !> \end{itemize}

    use glimmer_log
    use glimmer_coordinates, only: coordsystem_allocate

    implicit none

    type(glide_global_type),intent(inout) :: model

    integer :: ewn,nsn,upn

    ! for simplicity, copy these values...

    ewn = model%general%ewn
    nsn = model%general%nsn
    upn = model%general%upn
    
    ! horizontal coordinates

    allocate(model%general%x0(ewn-1))!; model%general%x0 = 0.d0  ! velocity grid
    allocate(model%general%y0(nsn-1))!; model%general%y0 = 0.d0
    allocate(model%general%x1(ewn))!; model%general%x1 = 0.d0    ! ice grid (for scalars)
    allocate(model%general%y1(nsn))!; model%general%y1 = 0.d0

    ! vertical sigma coordinates
    ! If we already have sigma, don't reallocate

    if (associated(model%numerics%sigma)) then
       if (size(model%numerics%sigma) /= upn) then
          call write_log('Wrong number of sigma levels given',GM_FATAL)
       end if
    else
       allocate(model%numerics%sigma(upn))
    endif

    allocate(model%numerics%stagsigma(upn-1))
    allocate(model%numerics%stagwbndsigma(0:upn))  !MJH added (0:upn) as separate variable

    ! temperature arrays

    !NOTE: In the glide dycore (whichdycore = DYCORE_GLIDE), the temperature and 
    !       flow factor live on the unstaggered vertical grid, and extra rows and columns 
    !       (with indices 0:ewn+1, 0:nsn+1) are needed.
    !      In the Glissade dycore, the temperature and flow factor live on
    !       the staggered vertical grid, with temp and flwa defined at the
    !       center of each layer k = 1:upn-1.  The temperature (but not flwa)
    !       is defined at the upper surface (k = 0) and lower surface (k = upn).

    if (model%options%whichdycore == DYCORE_GLIDE) then
       allocate(model%temper%temp(upn,0:ewn+1,0:nsn+1))
       call coordsystem_allocate(model%general%ice_grid, upn, model%temper%flwa)
       call coordsystem_allocate(model%general%ice_grid, upn, model%temper%dissip)
    else    ! glissade dycore
       allocate(model%temper%temp(0:upn,1:ewn,1:nsn))
       ! tempunstag has the same horizontal grid as the glissade temp, but a
       ! vertical axis like the glide temp
       allocate(model%temper%tempunstag(upn,1:ewn,1:nsn))
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%temper%flwa)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%temper%dissip)
    endif

    ! MJH - Set temp and flwa to physically unrealistic values so we can tell later if 
    !       arrays were initialized correctly
    model%temper%temp(:,:,:) = unphys_val  ! large negative number
    model%temper%flwa(:,:,:) = unphys_val
    model%temper%dissip(:,:,:) = 0.d0
    if (associated(model%temper%tempunstag)) &
         model%temper%tempunstag(:,:,:) = unphys_val

    call coordsystem_allocate(model%general%ice_grid,  model%temper%bheatflx)
    call coordsystem_allocate(model%general%ice_grid,  model%temper%bwat)
    call coordsystem_allocate(model%general%velo_grid, model%temper%stagbwat)
    call coordsystem_allocate(model%general%ice_grid,  model%temper%bpmp)
    call coordsystem_allocate(model%general%velo_grid, model%temper%stagbpmp)
    call coordsystem_allocate(model%general%ice_grid,  model%temper%btemp)
    call coordsystem_allocate(model%general%velo_grid, model%temper%stagbtemp)
    call coordsystem_allocate(model%general%ice_grid,  model%temper%ucondflx)

    if (model%options%whichdycore == DYCORE_GLIDE) then   ! glide only
       call coordsystem_allocate(model%general%ice_grid, model%temper%bwatflx)
    else   ! glissade only
       call coordsystem_allocate(model%general%ice_grid, model%temper%bfricflx)
       call coordsystem_allocate(model%general%ice_grid, model%temper%lcondflx)
       call coordsystem_allocate(model%general%ice_grid, model%temper%dissipcol)
       ! water fraction and enthalpy live at the midpoint of each layer (with temp and flwa)
       ! enthalpy (like temp) is defined at the upper and lower surfaces as well
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%temper%waterfrac)
       allocate(model%temper%enthalpy(0:upn,1:ewn,1:nsn))
       model%temper%enthalpy(:,:,:) = 0.d0
    endif

    ! velocity arrays

    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%uvel)
    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%vvel)
    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%velnorm)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%uflx)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%vflx)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%bed_softness)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%btrc)
    call coordsystem_allocate(model%general%velo_grid, 2, model%velocity%btraction)
    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%resid_u)
    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%resid_v)
    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%rhs_u)
    call coordsystem_allocate(model%general%velo_grid, upn, model%velocity%rhs_v)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%uvel_2d)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%vvel_2d)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%ubas)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%vbas)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%uvel_mean)
    call coordsystem_allocate(model%general%velo_grid, model%velocity%vvel_mean)
    call coordsystem_allocate(model%general%ice_grid,  upn, model%velocity%wvel)

    ! The following are on the extended staggered grid, which is the same size as the ice grid.
    call coordsystem_allocate(model%general%ice_grid,  upn, model%velocity%uvel_extend)
    call coordsystem_allocate(model%general%ice_grid,  upn, model%velocity%vvel_extend)
    call coordsystem_allocate(model%general%ice_grid,  model%velocity%uvel_2d_extend)
    call coordsystem_allocate(model%general%ice_grid,  model%velocity%vvel_2d_extend)

    if (model%options%whichdycore == DYCORE_GLIDE) then
       call coordsystem_allocate(model%general%ice_grid,  upn, model%velocity%wgrd)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%diffu)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%diffu_x)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%diffu_y)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%total_diffu)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%tau_x)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%tau_y)
    else   ! glissade dycore
       call coordsystem_allocate(model%general%velo_grid, model%velocity%beta)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%beta_internal)
       call coordsystem_allocate(model%general%ice_grid, model%velocity%unstagbeta)
       ! Set beta and unstagbeta to physically unrealistic values so we can tell later 
       ! if these fields were read correctly from an input file
       model%velocity%beta(:,:) = unphys_val   ! large negative number
       model%velocity%unstagbeta(:,:) = unphys_val

       call coordsystem_allocate(model%general%velo_grid, model%velocity%kinbcmask)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%umask_no_penetration)
       call coordsystem_allocate(model%general%velo_grid, model%velocity%vmask_no_penetration)
       call coordsystem_allocate(model%general%ice_grid, model%velocity%divu)
    endif

    ! higher-order stress arrays

    if (model%options%whichdycore /= DYCORE_GLIDE) then   ! glissade dycore
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%efvs)
       call coordsystem_allocate(model%general%ice_grid, model%stress%efvs_vertavg)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%tau%scalar) 
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%tau%xz)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%tau%yz)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%tau%xx)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%tau%yy)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%stress%tau%xy)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%velocity%strain_rate%scalar) 
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%velocity%strain_rate%xz)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%velocity%strain_rate%yz)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%velocity%strain_rate%xx)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%velocity%strain_rate%yy)
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%velocity%strain_rate%xy)
       call coordsystem_allocate(model%general%velo_grid, model%stress%btractx)
       call coordsystem_allocate(model%general%velo_grid, model%stress%btracty)
       call coordsystem_allocate(model%general%velo_grid, model%stress%btract)
       call coordsystem_allocate(model%general%ice_grid, model%stress%btractx_extend)
       call coordsystem_allocate(model%general%ice_grid, model%stress%btracty_extend)
       call coordsystem_allocate(model%general%velo_grid, model%stress%taudx)
       call coordsystem_allocate(model%general%velo_grid, model%stress%taudy)
    endif

    ! geometry arrays
    call coordsystem_allocate(model%general%ice_grid, model%geometry%thck)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%usrf)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%lsrf)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%topg)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%thck_obs)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%usrf_obs)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%topg_obs)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%thkmask)
    call coordsystem_allocate(model%general%velo_grid, model%geometry%stagmask)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%cell_area)

    call coordsystem_allocate(model%general%velo_grid, model%geomderv%stagthck)
    call coordsystem_allocate(model%general%velo_grid, model%geomderv%dthckdew)
    call coordsystem_allocate(model%general%velo_grid, model%geomderv%dthckdns)
    call coordsystem_allocate(model%general%velo_grid, model%geomderv%dusrfdew)
    call coordsystem_allocate(model%general%velo_grid, model%geomderv%dusrfdns)

    call coordsystem_allocate(model%general%ice_grid, model%geometry%sfc_mbal_flux)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%sfc_mbal_flux_tavg)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%basal_mbal_flux)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%basal_mbal_flux_tavg)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%calving_flux)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%calving_flux_tavg)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%gl_flux_east)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%gl_flux_north)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%gl_flux)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%gl_flux_tavg)

    call coordsystem_allocate(model%general%ice_grid, model%geometry%ice_mask)
    call coordsystem_allocate(model%general%velo_grid, model%geometry%ice_mask_stag)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%floating_mask)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%grounded_mask)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%ice_sheet_mask)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%ice_cap_mask)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%lower_cell_loc)
    call coordsystem_allocate(model%general%ice_grid, model%geometry%lower_cell_temp)

    if (model%options%whichdycore == DYCORE_GLIDE) then
       call coordsystem_allocate(model%general%ice_grid, model%geometry%thck_index)
       call coordsystem_allocate(model%general%ice_grid,  model%geomderv%dthckdtm)
       call coordsystem_allocate(model%general%ice_grid,  model%geomderv%dusrfdtm)
       allocate(model%thckwk%olds(ewn,nsn,model%thckwk%nwhich))
       model%thckwk%olds = 0.0d0
       call coordsystem_allocate(model%general%ice_grid, model%thckwk%oldthck)
       call coordsystem_allocate(model%general%ice_grid, model%thckwk%oldthck2)
    else   ! glissade dycore
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%geometry%ice_age)
       call coordsystem_allocate(model%general%ice_grid,  model%geometry%thck_old)
       call coordsystem_allocate(model%general%ice_grid,  model%geometry%dthck_dt)
       call coordsystem_allocate(model%general%ice_grid,  model%geometry%dthck_dt_tavg)
       call coordsystem_allocate(model%general%ice_grid,  model%geometry%f_flotation)
       call coordsystem_allocate(model%general%velo_grid, model%geometry%f_ground)
       call coordsystem_allocate(model%general%ice_grid,  model%geometry%f_ground_cell)
       call coordsystem_allocate(model%general%velo_grid, model%geomderv%dlsrfdew)
       call coordsystem_allocate(model%general%velo_grid, model%geomderv%dlsrfdns)
       call coordsystem_allocate(model%general%velo_grid, model%geomderv%staglsrf)
       call coordsystem_allocate(model%general%velo_grid, model%geomderv%stagusrf)
       call coordsystem_allocate(model%general%velo_grid, model%geomderv%stagtopg)
       !Note: model%geometry%tracers and related arrays are allocated later, in glissade_transport_setup

       ! Basal Physics
       !WHL - Since the number of basal BC options is proliferating, simplify the logic by allocating the following arrays
       !      whenever running glissade
!!       if ( (model%options%which_ho_babc == HO_BABC_POWERLAW) .or. &
!!            (model%options%which_ho_babc == HO_BABC_COULOMB_FRICTION) .or. &
!!            (model%options%which_ho_babc == HO_BABC_COULOMB_POWERLAW_SCHOOF) .or. &
!!            (model%options%whichbwat == BWATER_OCEAN_PENETRATION)     ) then
       call coordsystem_allocate(model%general%velo_grid, model%basal_physics%bpmp_mask)
       call coordsystem_allocate(model%general%ice_grid, model%basal_physics%effecpress)
       call coordsystem_allocate(model%general%velo_grid, model%basal_physics%effecpress_stag)
       call coordsystem_allocate(model%general%velo_grid, model%basal_physics%tau_c)
       call coordsystem_allocate(model%general%ice_grid, model%basal_physics%C_space_factor)
       call coordsystem_allocate(model%general%velo_grid, model%basal_physics%C_space_factor_stag)
       call coordsystem_allocate(model%general%velo_grid, model%basal_physics%mintauf)
!!       endif
    endif  ! glissade

    ! bmlt arrays
    call coordsystem_allocate(model%general%ice_grid,  model%basal_melt%bmlt)
    call coordsystem_allocate(model%general%ice_grid,  model%basal_melt%bmlt_applied)
    call coordsystem_allocate(model%general%ice_grid,  model%basal_melt%bmlt_applied_tavg)

    !WHL - debug
    call coordsystem_allocate(model%general%ice_grid,  model%basal_melt%bmlt_applied_old)
    call coordsystem_allocate(model%general%ice_grid,  model%basal_melt%bmlt_applied_diff)

    if (model%options%whichdycore == DYCORE_GLISSADE) then
       call coordsystem_allocate(model%general%ice_grid,  model%basal_melt%bmlt_ground)
       call coordsystem_allocate(model%general%ice_grid, model%basal_melt%bmlt_float)
       call coordsystem_allocate(model%general%ice_grid, model%basal_melt%bmlt_float_anomaly)
       call coordsystem_allocate(model%general%ice_grid, model%basal_melt%warm_ocean_mask)
       if (model%options%whichbmlt_float == BMLT_FLOAT_EXTERNAL) then
          call coordsystem_allocate(model%general%ice_grid, model%basal_melt%bmlt_float_external)
       endif
       if (model%options%whichbmlt_float == BMLT_FLOAT_MISOMIP) then
          call coordsystem_allocate(model%general%ice_grid, model%plume%T_basal)
          call coordsystem_allocate(model%general%ice_grid, model%plume%S_basal)
          call coordsystem_allocate(model%general%ice_grid, model%plume%u_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%v_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%u_plume_Cgrid)
          call coordsystem_allocate(model%general%ice_grid, model%plume%v_plume_Cgrid)
          call coordsystem_allocate(model%general%ice_grid, model%plume%D_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%ustar_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%drho_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%T_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%S_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%entrainment)
          call coordsystem_allocate(model%general%ice_grid, model%plume%detrainment)
          call coordsystem_allocate(model%general%ice_grid, model%plume%divDu_plume)
          call coordsystem_allocate(model%general%ice_grid, model%plume%T_ambient)
          call coordsystem_allocate(model%general%ice_grid, model%plume%S_ambient)
       endif
    endif  ! Glissade

    ! inversion arrays (Glissade only)
    if (model%options%which_ho_inversion == HO_INVERSION_COMPUTE .or. &
        model%options%which_ho_inversion == HO_INVERSION_PRESCRIBE) then
       call coordsystem_allocate(model%general%ice_grid, model%inversion%bmlt_float_inversion)
       call coordsystem_allocate(model%general%ice_grid, model%inversion%bmlt_float_prescribed)
       call coordsystem_allocate(model%general%ice_grid, model%inversion%powerlaw_c_inversion)
       call coordsystem_allocate(model%general%ice_grid, model%inversion%powerlaw_c_prescribed)
       call coordsystem_allocate(model%general%ice_grid, model%inversion%usrf_inversion)
       call coordsystem_allocate(model%general%ice_grid, model%inversion%dthck_dt_inversion)
    else
       ! Allocate powerlaw_c_inversion with size 1, since it is passed into subroutine calcbeta
       !  by the velocity solver.
       allocate(model%inversion%powerlaw_c_inversion(1,1))
    endif

    ! climate arrays
    call coordsystem_allocate(model%general%ice_grid, model%climate%acab)
    call coordsystem_allocate(model%general%ice_grid, model%climate%acab_tavg)
    call coordsystem_allocate(model%general%ice_grid, model%climate%acab_anomaly)
    call coordsystem_allocate(model%general%ice_grid, model%climate%acab_corrected)
    call coordsystem_allocate(model%general%ice_grid, model%climate%acab_applied)
    call coordsystem_allocate(model%general%ice_grid, model%climate%acab_applied_tavg)
    call coordsystem_allocate(model%general%ice_grid, model%climate%artm)
    call coordsystem_allocate(model%general%ice_grid, model%climate%smb)
    call coordsystem_allocate(model%general%ice_grid, model%climate%no_advance_mask)
    call coordsystem_allocate(model%general%ice_grid, model%climate%overwrite_acab_mask)

    ! calving arrays
    call coordsystem_allocate(model%general%ice_grid, model%calving%calving_thck)
    call coordsystem_allocate(model%general%ice_grid, model%calving%calving_rate)
    call coordsystem_allocate(model%general%ice_grid, model%calving%calving_rate_tavg)
    call coordsystem_allocate(model%general%ice_grid, model%calving%calving_mask)
    call coordsystem_allocate(model%general%ice_grid, model%calving%lateral_rate)
    call coordsystem_allocate(model%general%ice_grid, model%calving%tau_eigen1)
    call coordsystem_allocate(model%general%ice_grid, model%calving%tau_eigen2)
    call coordsystem_allocate(model%general%ice_grid, model%calving%tau_eff)
    if (model%options%whichcalving == CALVING_DAMAGE) then
       call coordsystem_allocate(model%general%ice_grid, upn-1, model%calving%damage)
    else
       ! allocate with size 1, since they need to be allocated to be passed to calving subroutine
       allocate(model%calving%damage(1,1,1))
    endif

    ! matrix solver arrays

    allocate (model%solver_data%rhsd(ewn*nsn))
    allocate (model%solver_data%answ(ewn*nsn))

    call new_sparse_matrix(ewn*nsn, 5*ewn*nsn, model%solver_data%matrix)

    !TODO - In model%lithot%temp, put the vertical index 3rd as in model%temper%temp?

    ! lithosphere arrays

    if (model%options%gthf == GTHF_COMPUTE) then
       allocate(model%lithot%temp(1:ewn,1:nsn,model%lithot%nlayer)); model%lithot%temp = 0.d0
       call coordsystem_allocate(model%general%ice_grid, model%lithot%mask)
    endif

    ! isostasy arrays

    ! Note: relx needs to be allocated always
    call coordsystem_allocate(model%general%ice_grid, model%isostasy%relx)
    if (model%options%isostasy == ISOSTASY_COMPUTE) then
       call coordsystem_allocate(model%general%ice_grid, model%isostasy%load)
       call coordsystem_allocate(model%general%ice_grid, model%isostasy%load_factors)
    endif

    ! The remaining arrays are not currently used
    ! phaml arrays
!!    call coordsystem_allocate(model%general%ice_grid, model%phaml%init_phaml)
!!    call coordsystem_allocate(model%general%ice_grid, model%phaml%rs_phaml)
!!    call coordsystem_allocate(model%general%ice_grid, model%phaml%uphaml)

    ! grounding line arrays (not currently supported)

!!    if (model%options%whichdycore /= DYCORE_GLIDE) then   ! glissade dycore
!!       allocate (model%ground%gl_ew(ewn-1,nsn))
!!       allocate (model%ground%gl_ns(ewn,nsn-1))
!!       allocate (model%ground%gline_flux(ewn,nsn)) 
!!    endif

    ! basal process arrays
    ! not currently supported

    if (model%options%whichdycore /= DYCORE_GLIDE) then   ! glissade dycore
!!       call coordsystem_allocate(model%general%ice_grid, model%basalproc%Hwater)
!!       allocate(model%basalproc%u (ewn-1,nsn-1,model%basalproc%tnodes)); model%basalproc%u=41.0d3
!!       allocate(model%basalproc%etill (ewn-1,nsn-1,model%basalproc%tnodes));model%basalproc%etill=0.5d0
    endif

  end subroutine glide_allocarr


  subroutine glide_deallocarr(model)

    !> deallocate model arrays
    !TODO - Check that all arrays allocated above are deallocated here.

    implicit none
    type(glide_global_type),intent(inout) :: model

    ! horizontal coordinates

    if (associated(model%general%x0)) &
        deallocate(model%general%x0) 
    if (associated(model%general%y0)) &
        deallocate(model%general%y0) 
    if (associated(model%general%x1)) &
        deallocate(model%general%x1) 
    if (associated(model%general%y1)) &
        deallocate(model%general%y1) 

    ! vertical sigma coordinates

    if (associated(model%numerics%sigma)) &
        deallocate(model%numerics%sigma)
    if (associated(model%numerics%stagsigma)) &
        deallocate(model%numerics%stagsigma)
    if (associated(model%numerics%stagwbndsigma)) &
        deallocate(model%numerics%stagwbndsigma)

    ! temperature arrays

    if (associated(model%temper%temp)) &
        deallocate(model%temper%temp)
    if (associated(model%temper%tempunstag)) &
        deallocate(model%temper%tempunstag)
    if (associated(model%temper%bheatflx)) &
        deallocate(model%temper%bheatflx)
    if (associated(model%temper%bwat)) &
        deallocate(model%temper%bwat)
    if (associated(model%temper%bwatflx)) &
        deallocate(model%temper%bwatflx)
    if (associated(model%temper%stagbwat)) &
        deallocate(model%temper%stagbwat)
    if (associated(model%temper%bpmp)) &
        deallocate(model%temper%bpmp)
    if (associated(model%temper%stagbpmp)) &
        deallocate(model%temper%stagbpmp)
    if (associated(model%temper%btemp)) &
        deallocate(model%temper%btemp)
    if (associated(model%temper%stagbtemp)) &
        deallocate(model%temper%stagbtemp)
    if (associated(model%temper%bfricflx)) &
        deallocate(model%temper%bfricflx)
    if (associated(model%temper%ucondflx)) &
        deallocate(model%temper%ucondflx)
    if (associated(model%temper%lcondflx)) &
        deallocate(model%temper%lcondflx)
    if (associated(model%temper%dissipcol)) &
        deallocate(model%temper%dissipcol)
    if (associated(model%temper%waterfrac)) &
        deallocate(model%temper%waterfrac)
    if (associated(model%temper%enthalpy)) &
        deallocate(model%temper%enthalpy)
    if (associated(model%temper%flwa)) &
        deallocate(model%temper%flwa)
    if (associated(model%temper%dissip)) &
        deallocate(model%temper%dissip)

    ! velocity arrays

    if (associated(model%velocity%uvel)) &
        deallocate(model%velocity%uvel)
    if (associated(model%velocity%vvel)) &
        deallocate(model%velocity%vvel)
    if (associated(model%velocity%uvel_2d)) &
        deallocate(model%velocity%uvel_2d)
    if (associated(model%velocity%vvel_2d)) &
        deallocate(model%velocity%vvel_2d)
    if (associated(model%velocity%velnorm)) &
        deallocate(model%velocity%velnorm)
    if (associated(model%velocity%wvel)) &
        deallocate(model%velocity%wvel)
    if (associated(model%velocity%uflx)) &
        deallocate(model%velocity%uflx)
    if (associated(model%velocity%vflx)) &
        deallocate(model%velocity%vflx)
    if (associated(model%velocity%bed_softness)) &
        deallocate(model%velocity%bed_softness)
    if (associated(model%velocity%btrc)) &
        deallocate(model%velocity%btrc)
    if (associated(model%velocity%btraction)) &
        deallocate(model%velocity%btraction)
    if (associated(model%velocity%resid_u)) &
        deallocate(model%velocity%resid_u)
    if (associated(model%velocity%resid_v)) &
        deallocate(model%velocity%resid_v)
    if (associated(model%velocity%rhs_u)) &
        deallocate(model%velocity%rhs_u)
    if (associated(model%velocity%rhs_v)) &
        deallocate(model%velocity%rhs_v)
    if (associated(model%velocity%uvel_extend)) &
        deallocate(model%velocity%uvel_extend)
    if (associated(model%velocity%vvel_extend)) &
        deallocate(model%velocity%vvel_extend)
    if (associated(model%velocity%uvel_2d_extend)) &
        deallocate(model%velocity%uvel_2d_extend)
    if (associated(model%velocity%vvel_2d_extend)) &
        deallocate(model%velocity%vvel_2d_extend)
    if (associated(model%velocity%ubas)) &
        deallocate(model%velocity%ubas)
    if (associated(model%velocity%vbas)) &
        deallocate(model%velocity%vbas)
    if (associated(model%velocity%uvel_mean)) &
        deallocate(model%velocity%uvel_mean)
    if (associated(model%velocity%vvel_mean)) &
        deallocate(model%velocity%vvel_mean)

    if (associated(model%velocity%wgrd)) &
        deallocate(model%velocity%wgrd)
    if (associated(model%velocity%diffu)) &
        deallocate(model%velocity%diffu)
    if (associated(model%velocity%diffu_x)) &
        deallocate(model%velocity%diffu_x)
    if (associated(model%velocity%diffu_y)) &
        deallocate(model%velocity%diffu_y)
    if (associated(model%velocity%total_diffu)) &
        deallocate(model%velocity%total_diffu)
    if (associated(model%velocity%tau_x)) &
        deallocate(model%velocity%tau_x)
    if (associated(model%velocity%tau_y)) &
        deallocate(model%velocity%tau_y)

!!    if (associated(model%velocity%velmask)) &   ! no longer used
!!       deallocate(model%velocity%velmask) 

    if (associated(model%velocity%beta)) &
        deallocate(model%velocity%beta)
    if (associated(model%velocity%unstagbeta)) &
        deallocate(model%velocity%unstagbeta)
    if (associated(model%velocity%beta_internal)) &
        deallocate(model%velocity%beta_internal)
    if (associated(model%velocity%kinbcmask)) &
        deallocate(model%velocity%kinbcmask)
    if (associated(model%velocity%umask_no_penetration)) &
        deallocate(model%velocity%umask_no_penetration)
    if (associated(model%velocity%vmask_no_penetration)) &
        deallocate(model%velocity%vmask_no_penetration)
    if (associated(model%velocity%divu)) &
        deallocate(model%velocity%divu)

    ! higher-order stress arrays

    if (associated(model%stress%efvs)) &
        deallocate(model%stress%efvs)
    if (associated(model%stress%efvs_vertavg)) &
        deallocate(model%stress%efvs_vertavg)
    if (associated(model%stress%tau%scalar)) &
        deallocate(model%stress%tau%scalar)
    if (associated(model%stress%tau%xz)) &
        deallocate(model%stress%tau%xz)
    if (associated(model%stress%tau%yz)) &
        deallocate(model%stress%tau%yz)
    if (associated(model%stress%tau%xx)) &
        deallocate(model%stress%tau%xx)
    if (associated(model%stress%tau%yy)) &
        deallocate(model%stress%tau%yy)
    if (associated(model%stress%tau%xy)) &
        deallocate(model%stress%tau%xy)
    if (associated(model%velocity%strain_rate%scalar)) &
        deallocate(model%velocity%strain_rate%scalar)
    if (associated(model%velocity%strain_rate%xz)) &
        deallocate(model%velocity%strain_rate%xz)
    if (associated(model%velocity%strain_rate%yz)) &
        deallocate(model%velocity%strain_rate%yz)
    if (associated(model%velocity%strain_rate%xx)) &
        deallocate(model%velocity%strain_rate%xx)
    if (associated(model%velocity%strain_rate%yy)) &
        deallocate(model%velocity%strain_rate%yy)
    if (associated(model%velocity%strain_rate%xy)) &
        deallocate(model%velocity%strain_rate%xy)
    if (associated(model%stress%btractx)) &
        deallocate(model%stress%btractx)
    if (associated(model%stress%btracty)) &
        deallocate(model%stress%btracty)
    if (associated(model%stress%btract)) &
        deallocate(model%stress%btract)
    if (associated(model%stress%btractx_extend)) &
        deallocate(model%stress%btractx_extend)
    if (associated(model%stress%btracty_extend)) &
        deallocate(model%stress%btracty_extend)
    if (associated(model%stress%taudx)) &
        deallocate(model%stress%taudx)
    if (associated(model%stress%taudy)) &
        deallocate(model%stress%taudy)

    ! basal physics arrays
    if (associated(model%basal_physics%bpmp_mask)) &
        deallocate(model%basal_physics%bpmp_mask)
    if (associated(model%basal_physics%effecpress)) &
        deallocate(model%basal_physics%effecpress)
    if (associated(model%basal_physics%effecpress_stag)) &
        deallocate(model%basal_physics%effecpress_stag)
    if (associated(model%basal_physics%tau_c)) &
        deallocate(model%basal_physics%tau_c)
    if (associated(model%basal_physics%C_space_factor)) &
        deallocate(model%basal_physics%C_space_factor)
    if (associated(model%basal_physics%C_space_factor_stag)) &
        deallocate(model%basal_physics%C_space_factor_stag)
    if (associated(model%basal_physics%mintauf)) &
       deallocate(model%basal_physics%mintauf)

    ! basal melt arrays

    if (associated(model%basal_melt%bmlt)) &
        deallocate(model%basal_melt%bmlt)
    if (associated(model%basal_melt%bmlt_applied)) &
        deallocate(model%basal_melt%bmlt_applied)
    if (associated(model%basal_melt%bmlt_applied_tavg)) &
        deallocate(model%basal_melt%bmlt_applied_tavg)
    if (associated(model%basal_melt%bmlt_ground)) &
        deallocate(model%basal_melt%bmlt_ground)
    if (associated(model%basal_melt%bmlt_float)) &
        deallocate(model%basal_melt%bmlt_float)
    if (associated(model%basal_melt%bmlt_float_external)) &
        deallocate(model%basal_melt%bmlt_float_external)
    if (associated(model%basal_melt%bmlt_float_anomaly)) &
        deallocate(model%basal_melt%bmlt_float_anomaly)
    if (associated(model%basal_melt%warm_ocean_mask)) &
        deallocate(model%basal_melt%warm_ocean_mask)
    if (associated(model%basal_melt%bmlt_applied_old)) &
        deallocate(model%basal_melt%bmlt_applied_old)
    if (associated(model%basal_melt%bmlt_applied_diff)) &
        deallocate(model%basal_melt%bmlt_applied_diff)

    ! inversion arrays
    if (associated(model%inversion%bmlt_float_inversion)) &
        deallocate(model%inversion%bmlt_float_inversion)
    if (associated(model%inversion%bmlt_float_prescribed)) &
        deallocate(model%inversion%bmlt_float_prescribed)
    if (associated(model%inversion%powerlaw_c_inversion)) &
        deallocate(model%inversion%powerlaw_c_inversion)
    if (associated(model%inversion%powerlaw_c_prescribed)) &
        deallocate(model%inversion%powerlaw_c_prescribed)
    if (associated(model%inversion%usrf_inversion)) &
        deallocate(model%inversion%usrf_inversion)
    if (associated(model%inversion%dthck_dt_inversion)) &
        deallocate(model%inversion%dthck_dt_inversion)

    ! plume arrays
    if (associated(model%plume%T_basal)) &
        deallocate(model%plume%T_basal)
    if (associated(model%plume%S_basal)) &
        deallocate(model%plume%S_basal)
    if (associated(model%plume%u_plume)) &
        deallocate(model%plume%u_plume)
    if (associated(model%plume%v_plume)) &
        deallocate(model%plume%v_plume)
    if (associated(model%plume%u_plume_Cgrid)) &
        deallocate(model%plume%u_plume_Cgrid)
    if (associated(model%plume%v_plume_Cgrid)) &
        deallocate(model%plume%v_plume_Cgrid)
    if (associated(model%plume%D_plume)) &
        deallocate(model%plume%D_plume)
    if (associated(model%plume%ustar_plume)) &
        deallocate(model%plume%ustar_plume)
    if (associated(model%plume%drho_plume)) &
        deallocate(model%plume%drho_plume)
    if (associated(model%plume%T_plume)) &
        deallocate(model%plume%T_plume)
    if (associated(model%plume%S_plume)) &
        deallocate(model%plume%S_plume)
    if (associated(model%plume%entrainment)) &
        deallocate(model%plume%entrainment)
    if (associated(model%plume%detrainment)) &
        deallocate(model%plume%detrainment)
    if (associated(model%plume%divDu_plume)) &
        deallocate(model%plume%divDu_plume)
    if (associated(model%plume%T_ambient)) &
        deallocate(model%plume%T_ambient)
    if (associated(model%plume%S_ambient)) &
        deallocate(model%plume%S_ambient)

    ! geometry arrays

    if (associated(model%geometry%thck)) &
        deallocate(model%geometry%thck)
    if (associated(model%geometry%usrf)) &
        deallocate(model%geometry%usrf)
    if (associated(model%geometry%lsrf)) &
        deallocate(model%geometry%lsrf)
    if (associated(model%geometry%topg)) &
        deallocate(model%geometry%topg)
    if (associated(model%geometry%thck_obs)) &
        deallocate(model%geometry%thck_obs)
    if (associated(model%geometry%usrf_obs)) &
        deallocate(model%geometry%usrf_obs)
    if (associated(model%geometry%topg_obs)) &
        deallocate(model%geometry%topg_obs)
    if (associated(model%geometry%thkmask)) &
        deallocate(model%geometry%thkmask)
    if (associated(model%geometry%stagmask)) &
        deallocate(model%geometry%stagmask)
    if (associated(model%geomderv%stagthck)) &
        deallocate(model%geomderv%stagthck)
    if (associated(model%geomderv%dthckdew)) &
        deallocate(model%geomderv%dthckdew)
    if (associated(model%geomderv%dthckdns)) &
        deallocate(model%geomderv%dthckdns)
    if (associated(model%geomderv%dusrfdew)) &
        deallocate(model%geomderv%dusrfdew)
    if (associated(model%geomderv%dusrfdns)) &
        deallocate(model%geomderv%dusrfdns)

    if (associated(model%geometry%cell_area)) &
        deallocate(model%geometry%cell_area)

    if (associated(model%geometry%sfc_mbal_flux)) &
        deallocate(model%geometry%sfc_mbal_flux)
    if (associated(model%geometry%sfc_mbal_flux_tavg)) &
        deallocate(model%geometry%sfc_mbal_flux_tavg)
    if (associated(model%geometry%basal_mbal_flux)) &
        deallocate(model%geometry%basal_mbal_flux)
    if (associated(model%geometry%basal_mbal_flux_tavg)) &
        deallocate(model%geometry%basal_mbal_flux_tavg)
    if (associated(model%geometry%calving_flux)) &
        deallocate(model%geometry%calving_flux)
    if (associated(model%geometry%calving_flux_tavg)) &
        deallocate(model%geometry%calving_flux_tavg)
    if (associated(model%geometry%gl_flux_east)) &
        deallocate(model%geometry%gl_flux_east)
    if (associated(model%geometry%gl_flux_north)) &
        deallocate(model%geometry%gl_flux_north)
    if (associated(model%geometry%gl_flux)) &
        deallocate(model%geometry%gl_flux)
    if (associated(model%geometry%gl_flux_tavg)) &
        deallocate(model%geometry%gl_flux_tavg)

    if (associated(model%geometry%ice_mask)) &
       deallocate(model%geometry%ice_mask)
    if (associated(model%geometry%ice_mask_stag)) &
       deallocate(model%geometry%ice_mask_stag)
    if (associated(model%geometry%floating_mask)) &
       deallocate(model%geometry%floating_mask)
    if (associated(model%geometry%grounded_mask)) &
       deallocate(model%geometry%grounded_mask)
    if (associated(model%geometry%ice_sheet_mask)) &
       deallocate(model%geometry%ice_sheet_mask)
    if (associated(model%geometry%ice_cap_mask)) &
       deallocate(model%geometry%ice_cap_mask)
    if (associated(model%geometry%lower_cell_loc)) &
       deallocate(model%geometry%lower_cell_loc)
    if (associated(model%geometry%lower_cell_temp)) &
       deallocate(model%geometry%lower_cell_temp)

    if (associated(model%geometry%thck_index)) &
        deallocate(model%geometry%thck_index)
    if (associated(model%geomderv%dthckdtm)) &
        deallocate(model%geomderv%dthckdtm)
    if (associated(model%geomderv%dusrfdtm)) &
        deallocate(model%geomderv%dusrfdtm)
    if (associated(model%thckwk%olds)) &
        deallocate(model%thckwk%olds)
    if (associated(model%thckwk%oldthck)) &
        deallocate(model%thckwk%oldthck)
    if (associated(model%thckwk%oldthck2)) &
        deallocate(model%thckwk%oldthck2)

    if (associated(model%geometry%ice_age)) &
        deallocate(model%geometry%ice_age)
    if (associated(model%geometry%thck_old)) &
        deallocate(model%geometry%thck_old)
    if (associated(model%geometry%dthck_dt)) &
        deallocate(model%geometry%dthck_dt)
    if (associated(model%geometry%dthck_dt_tavg)) &
        deallocate(model%geometry%dthck_dt_tavg)
    if (associated(model%geometry%tracers)) &
        deallocate(model%geometry%tracers)
    if (associated(model%geometry%f_flotation)) &
        deallocate(model%geometry%f_flotation)
    if (associated(model%geometry%f_ground)) &
        deallocate(model%geometry%f_ground)
    if (associated(model%geometry%f_ground_cell)) &
        deallocate(model%geometry%f_ground_cell)
    if (associated(model%geomderv%dlsrfdew)) &
        deallocate(model%geomderv%dlsrfdew)
    if (associated(model%geomderv%dlsrfdns)) &
        deallocate(model%geomderv%dlsrfdns)
    if (associated(model%geomderv%staglsrf)) &
        deallocate(model%geomderv%staglsrf)
    if (associated(model%geomderv%stagusrf)) &
        deallocate(model%geomderv%stagusrf)
    if (associated(model%geomderv%stagtopg)) &
        deallocate(model%geomderv%stagtopg)

    ! climate arrays

    if (associated(model%climate%acab)) &
        deallocate(model%climate%acab)
    if (associated(model%climate%acab_tavg)) &
        deallocate(model%climate%acab_tavg)
    if (associated(model%climate%acab_anomaly)) &
        deallocate(model%climate%acab_anomaly)
    if (associated(model%climate%acab_corrected)) &
        deallocate(model%climate%acab_corrected)
    if (associated(model%climate%acab_applied)) &
        deallocate(model%climate%acab_applied)
    if (associated(model%climate%acab_applied_tavg)) &
        deallocate(model%climate%acab_applied_tavg)
    if (associated(model%climate%smb)) &
        deallocate(model%climate%smb)
    if (associated(model%climate%artm)) &
        deallocate(model%climate%artm)
    if (associated(model%climate%no_advance_mask)) &
        deallocate(model%climate%no_advance_mask)
    if (associated(model%climate%overwrite_acab_mask)) &
        deallocate(model%climate%overwrite_acab_mask)

    ! calving arrays
    if (associated(model%calving%calving_thck)) &
        deallocate(model%calving%calving_thck)
    if (associated(model%calving%calving_rate)) &
        deallocate(model%calving%calving_rate)
    if (associated(model%calving%calving_rate_tavg)) &
        deallocate(model%calving%calving_rate_tavg)
    if (associated(model%calving%calving_mask)) &
        deallocate(model%calving%calving_mask)
    if (associated(model%calving%lateral_rate)) &
        deallocate(model%calving%lateral_rate)
    if (associated(model%calving%tau_eigen1)) &
        deallocate(model%calving%tau_eigen1)
    if (associated(model%calving%tau_eigen2)) &
        deallocate(model%calving%tau_eigen2)
    if (associated(model%calving%tau_eff)) &
        deallocate(model%calving%tau_eff)
    if (associated(model%calving%damage)) &
        deallocate(model%calving%damage)

    ! matrix solver arrays

    if (associated(model%solver_data%rhsd))  &  
        deallocate(model%solver_data%rhsd)
    if (associated(model%solver_data%answ))  &
        deallocate(model%solver_data%answ)

    ! lithosphere arrays

    if (associated(model%lithot%temp)) &
        deallocate(model%lithot%temp)
    if (associated(model%lithot%mask)) &
        deallocate(model%lithot%mask)

    ! isostasy arrays

    if (associated(model%isostasy%relx)) &
        deallocate(model%isostasy%relx)
    if (associated(model%isostasy%load)) &
        deallocate(model%isostasy%load)
    if (associated(model%isostasy%load_factors)) &
        deallocate(model%isostasy%load_factors)

    ! projection arrays
    if (associated(model%projection%stere)) then
       if (associated(model%projection%stere%area_factor)) then
          deallocate(model%projection%stere%area_factor)
       endif
    endif

    ! The remaining arrays are not currently used
    ! phaml arrays

!!    if (associated(model%phaml%init_phaml)) &
!!       deallocate(model%phaml%init_phaml)
!!    if (associated(model%phaml%rs_phaml)) &
!!       deallocate(model%phaml%rs_phaml)    
!!    if (associated(model%phaml%uphaml)) &
!!       deallocate(model%phaml%uphaml)

    ! grounding line arrays (not currently supported)

!!    if (associated(model%ground%gl_ns)) &
!!        deallocate(model%ground%gl_ns)
!!    if (associated(model%ground%gl_ew)) &
!!        deallocate(model%ground%gl_ew)
!!    if (associated(model%ground%gline_flux)) &
!!        deallocate(model%ground%gline_flux)

    ! basal process arrays
    ! not currently supported

!!    if (associated(model%basalproc%Hwater)) &
!!       deallocate(model%basalproc%Hwater)
!!    if (associated(model%basalproc%u)) &
!!       deallocate(model%basalproc%u)
!!    if (associated(model%basalproc%etill)) &
!!       deallocate(model%basalproc%etill)

  end subroutine glide_deallocarr


  ! some accessor functions
  function get_dew(model)
    !> return scaled x node spacing
    use glimmer_paramets, only : len0
    implicit none
    real(dp) :: get_dew
    type(glide_global_type) :: model

    get_dew = model%numerics%dew * len0
  end function get_dew

  function get_dns(model)
    !> return scaled y node spacing
    use glimmer_paramets, only : len0
    implicit none
    real(dp) :: get_dns
    type(glide_global_type) :: model

    get_dns = model%numerics%dns * len0
  end function get_dns

  function get_tstart(model)
    !> return start time
    implicit none
    real(dp) :: get_tstart
    type(glide_global_type) :: model
    
    get_tstart = model%numerics%tstart
  end function get_tstart

  function get_tend(model)
    !> return end time
    implicit none
    real(dp) :: get_tend
    type(glide_global_type) :: model
    
    get_tend = model%numerics%tend
  end function get_tend

  function get_tinc(model)
    !> return time increment
    implicit none
    real(dp) :: get_tinc
    type(glide_global_type) :: model
    
    get_tinc = model%numerics%tinc
  end function get_tinc

  function get_ewn(model)
    !> get number of nodes in x dir
    implicit none
    integer get_ewn
    type(glide_global_type) :: model

    get_ewn = model%general%ewn
  end function get_ewn

  function get_nsn(model)
    !> get number of nodes in y dir
    implicit none
    integer get_nsn
    type(glide_global_type) :: model

    get_nsn = model%general%nsn
  end function get_nsn
  
  subroutine set_time(model,time)
    !> Set the model time counter --- useful for
    !> fractional year output
    implicit none
    type(glide_global_type) :: model
    real(dp) :: time

    model%numerics%time = time
  end subroutine set_time

end module glide_types
